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1.  Introduction 


This  report  documents  the  progress  during  the  first  year  of  the  Global  Ground 
Tomography  (GGT)  contract,  a  research  effort  in  imaging  of  underground  structures  using 
low  frequency  electromagnetic  (EM)  signals.  This  effort  employs  theoretical  analysis, 
simulation,  and  experiments  to  assess  the  performance  of  low  frequency  EM  imaging  and 
develop  an  understanding  of  the  basic  physical  relationships  that  would  affect  the  design 
of  an  operational  imaging  system.  This  report  provides  a  detailed  description  of  the  work 
performed  to  date  in  each  of  the  three  major  areas  included  in  the  contract:  signal 
processing  and  imaging,  electromagnetic  simulation,  and  experiments. 

In  the  signal  processing  area,  the  major  results  to  date  have  been  the  development  and 
implementation  of  software  for  computing  impedance  estimates  from  experimental  data; 
this  software  has  been  used  for  processing  all  of  the  data  acquired  under  this  contract. 
Work  in  the  imaging  area  has  concentrated  on  applying  APTI’s  proprietary  imaging 
techniques  to  both  simulated  and  actual  data,  and  comparing  the  results  to  those  obtained 
using  conventional  imaging  algorithms.  A  detailed  discussion  of  the  differences  between 
APTI’s  method  and  conventional  approaches  is  given  in  Section  2.2.  Only  preliminary 
results  of  the  comparisons  are  available,  and  some  of  these  are  shown  in  Section  3.2; 
more  detailed  direct  comparisons  using  both  simulated  and  actual  data  are  underway. 

Electromagnetic  simulations  of  targets  similar  to  those  employed  in  the  first  two  field 
experiments  have  been  computed  and  analyzed.  These  two-dimensional  simulations 
model  empty  (resistive)  tunnels  in  one-  or  two-layer  grounds.  Even  for  this  simplified 
case,  several  important  effects  have  been  observed,  these  are  discussed  in  Section  3.  The 
results  of  employing  a  widely  used  imaging  algorithm  to  simulated  data  are  also 
discussed  in  this  section.  Future  efforts  will  include  modeling  of  conductive  two- 
dimensional  targets  and  initial  examinations  of  realistic  three-dimensional  targets  such  as 
buried  facilities  and  bunkers. 
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Progress  in  the  experimental  area  included  the  design  and  integration  of  a  data  acquisition 
system  described  in  Section  4.1,  and  the  use  of  this  system  to  conduct  two  field 
experiments.  The  first  experiment  employed  a  local  VLF  source,  and  was  conducted  at 
the  Superconducting  Supercollider  site  near  Waxahachee,  Texas;  the  second  used  both  a 
local  transmitter  and  the  HAARP  ionospheric  heater  as  sources,  and  was  carried  out  at  a 
mine  near  Cantwell,  Alaska.  The  results  of  these  experiments  are  presented  and  discussed 
in  Section  4.2.  Although  not  all  of  the  data  has  yet  been  analyzed,  a  successful  result  has 
been  obtained  in  the  Alaskan  experiment,  locating  a  mine  tunnel  at  a  depth  of  29  meters. 

The  final  section  summarizes  the  results  and  outlines  plans  for  the  next  year  of  this  effort. 
Major  goals  include  elimination  of  the  backlog  of  unprocessed  experimental  data, 
improvements  in  the  speed  of  data  processing  and  imaging,  more  extensive  use  of 
simulations  in  detailed  studies  of  imaging  performance,  and  execution  and  analysis  of 
additional  field  experiments  on  a  wider  variety  of  targets. 
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2.  Signal  Processing  and  Imaging 


In  the  following  section,  an  outline  of  the  signal  processing  and  imaging  techniques 
applied  to  convert  measured  data  into  images  of  underground  conductivity  variations  is 
given.  This  process  begins  with  the  processing  of  the  measured  data  into  impedance 
estimates,  and  then  converts  the  impedance  estimates  from  several  locations  into  an 
image  of  the  underground  conductivity. 

2.1.  Signal  Processing 

The  first  step  in  processing  measured  data  is  to  convert  the  time  series  for  the  electric  and 
magnetic  fields  measured  at  a  given  location  into  an  estimate  of  the  surface  impedance  at 
that  location.  This  process  is  described  in  the  sections  below. 

2.1.1.  Impedance  Estimation 

From  measurements  of  electric  and  magnetic  fields  at  a  given  location,  we  would  like  to 
determine  the  surface  impedance  as  a  function  of  frequency  at  that  location.  The  basic 
equation  that  defines  the  surface  impedance  matrix  in  terms  of  the  electric  and  magnetic 
fields  is 


E»' 

£»_ 


All  of  the  entries  in  this  matrix  equation  are  functions  of  frequency,  so  all  of  the 
operations  we  will  describe  below  must  be  performed  for  each  frequency  of  interest. 
Often,  we  will  drop  the  explicit  dependence  on  frequency  to  make  the  equations  more 
compact. 


The  matrix  equation  that  defines  the  impedances  is  a  system  of  two  equations  in  four 
unknowns  that  cannot  be  solved  directly.  It  is  possible  to  produce  a  system  of  four 
equations  in  four  unknowns  by  multiplying  both  sides  of  the  equation  on  the  right  by  a 
two-element  row  vector,  and  then  computing  the  expected  value  of  the  result  (an 
operation  indicated  by  angle  brackets).  We  will  assume  that  the  statistical  expected  value 
is  the  same  as  the  value  obtained  by  averaging  the  quantities  over  an  infinitely  long  data 


3 


record.  Our  initial  analysis  will  then  indicate  the  results  that  would  be  obtained  with 
extremely  large  amounts  of  measured  data. 


If  we  choose  the  elements  of  the  row  vector  to  be  the  complex  conjugates  of  the  field 
values  Hx  and  Hy,  for  example,  we  obtain  the  equation 


£»' 

EM 


[h;(co)  h* (&>)])  = 


ZyM  ZyyH 


yy' 


nx(co) 

hM 


[h»  h»] 


Choosing  to  make  the  elements  of  the  row  vector  complex  conjugates  of  some  quantity 
makes  the  elements  of  the  resulting  matrices  auto-  or  cross-power  spectral  densities,  so 
that  the  elements  themselves  may  be  easily  interpreted.  For  example,  the  equation  above 
reduces  to 


$ExHx  (^)  ^ExHy  (^) 

x» 

z„W~ 

SfJxHx  (^) 

S HxHy  (^ ) 

EyHx (^)  SEyHy((0)- 

Zyx(C0) 

Z„(a)_ 

J^HyHx^t0) 

SfiyHy  (^) 

where  all  the  elements  to  be  computed  from  the  measured  data  are  power  spectral 
densities. 

For  the  resulting  equation  to  be  nontrivial  (i.e.  not  just  the  matrix  equation  0  =  0),  the 
elements  of  the  row  vector  should  be  correlated  statistically  with  the  field  measurements. 
This  requirement  is  clearly  satisfied  by  conjugates  of  the  measurements  themselves,  or  of 
any  similar  measurement  of  electric  or  magnetic  fields  at  nearby  locations.  The  two  most 
common  choices  are  to  employ  [#*  H*  j ,  as  in  the  example  above,  or  to  use  two  remote 

reference  signals  [/?*  7?*],  as  components  of  the  row  vector.  The  first  approach  is  called 

the  standard  method,  and  the  second,  the  remote  reference  method. 

Although  both  of  these  approaches  produce  similar  results  in  the  absence  of  noise,  when 
noise  is  present,  the  remote  reference  technique  has  significant  advantages.  To  see  why, 
let  us  examine  the  case  where  all  of  the  field  measurements  have  measurement  noise 
added  to  them,  for  example,  the  internal  noise  of  the  field  sensors.  This  changes  the 
standard  method  equation  above  to 
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£,  +  A^ 

Ey  +  NEy 


fc’+A^  n;  +  N'Hy] )  = 


l\_ 

N* 

1 _ 

_ 1 

/ 

Hx  +  NHx 

r 

_ZyX 

l  _  . 

\ 

Hy  + 

K  +  N'Hx 


If  we  assume  that  the  noise  in  each  measurement  is  statistically  uncorrelated  with  the 
noise  in  all  the  other  measurements,  the  standard  method  results  in  the  equation 


^  Ex  Hx 

^  ExHy 

Zxc 

V 

“c  +  W 

°HxHx  ~  YYHxHx 

S HxHy 

^EyHx 

^EyHy  _ 

_Zyx 

Zyy^ 

^HyHx 

S  HyHy  +  ^HyHy  _ 

where  WhxHx  and  WuyHy  are  the  power  spectral  densities  of  the  measurement  noise  added 
to  Hx  and  Hy.  This  is  clearly  a  different  equation  than  the  noise-free  version,  and  its 
solution  will  be  different.  To  see  how  the  addition  of  noise  affects  the  solution,  let  us 
examine  the  changes  in  one  component,  Z^.  The  noise-free  equation  has  the  solution 

_  ^ExHy^HxHx  ~  ^ExEx^HxHy 

W  ~  C  C  _  C  C  ’ 

°HxHx°HyHy  °HxHy°HyHx 


while  the  noisy  equation  has  the  different  solution 


z*y  5 


^ExHyi.^ 


+w 

HxHx  ^  yyHxHx 


)-s 


ExEx^HxHy 


HyHy  \u  HxHx  +  ^HxHx  )  ^HxHy^HyHx  +  ^  HxHx^HyHy  +  ^HxHx^HyHy 


To  simplify  these  two  equations,  assume  that  the  measurements  are  taken  over  a  one 
dimensional  vertically  stratified  region  with  no  horizontal  variation.  This  causes  several 
of  the  crosspower  spectral  densities  to  become  zero,  and  the  noise-free  equation  becomes 

$ ExHy 

■ xy  c  ’ 

°  HyHy 


while  the  noisy  equation  in  this  case  is 

2 _ ^ExHy _ 

^HyHy  +  HxHx^HyHy  +  ^HxHx^HyHy  )/  i^HxHx  +  ^HxHx  ) 

We  can  see  that  the  magnitude  of  the  denominator  of  the  noisy  equation  is  the  same  or 
larger,  because  all  of  the  extra  terms  are  autopower  spectral  densities,  which  must  be 
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nonnegative.  The  magnitude  of  the  answer  will  therefore  be  the  same  or  smaller.  The 
results  of  this  particular  case  also  hold  in  general:  the  use  of  the  standard  method  when 
measurement  noise  is  not  negligible  results  in  impedance  estimates  whose  magnitudes  are 
lower  than  their  true  values.  Note  that  this  is  true  even  in  the  limit  of  infinitely  long 
measurement  times;  this  type  of  error  cannot  be  averaged  out. 


In  contrast,  if  we  use  the  remote  reference  method,  the  equation  to  be  solved  is 


£»' 

£» 


IXV)  <(«)])= 


Z„(a>) 


ZyM) 


xy' 

Zyy(to) 


'«»■ 

hM 


K(®) 


which  reduces  in  the  noise-free  case  to 

^HxRa  (^)  SHxRb  (CO) 

^HyRa  ^ HyRb^\ 


$ ExRa (^)  $ ExRb (^) 

Zjco)  Zjco) 

_^EyRa^^)  ^  EyRb^\ 

Zyx{(o)  Zyy{co)_ 

When  noise  is  added  to  the  measurements,  the  equation  becomes 


Ex  +  NEx 

Ey  +  N  Ey 


iK  +  Nl  R*b+N;b]) 


\ 

X 

xl 

X  +  NHx 

r 

X 

\ 

_Hy  +  NHy_ 

[K  +  Ka  K  +  *rJj' 


which  reduces  to  the  same  equation  as  the  noise  free  case: 


S ExRa 

S ExRb 

^xx 

xl 

^HxRa 

S HxRb 

_^EyRa 

$ EyRb  _ 

N* 
_ 1 

1 _ 

^HyRb  _ 

The  perturbing  terms  are  absent  because  the  noise  in  each  of  the  remote  reference 
channels  Ra  and  Rt,  is  independent  of  the  noise  in  Hx  and  Hy.  Thus,  even  in  the  case  where 
noise  is  present,  we  can  obtain  accurate  estimates  of  the  surface  impedances  by  using  the 
remote  reference  method,  as  long  as  our  measurement  times  are  sufficiently  long. 

The  remote  reference  method  would  eliminate  the  bias  due  to  measurement  noise  even  if 
the  sensors  providing  the  remote  reference  signal  were  placed  at  the  same  location  as  the 
E  and  H  sensors.  Placing  the  remote  reference  sensors  at  a  distance  from  the  other  sensors 
brings  an  additional  advantage  by  allowing  the  technique  to  remove  bias  due  to  sources 
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other  than  the  internal  noise  of  the  sensor,  such  as  local  electromagnetic  interference  or 
vibration-induced  signals  in  magnetic  field  sensors.  As  long  as  the  remote  reference  is 
placed  so  that  any  signal  present  at  both  the  reference  location  and  the  measurement 
location  must  be  an  electromagnetic  wave  propagating  in  the  Earth-ionosphere 
waveguide,  the  results  obtained  with  this  technique  will  not  be  biased  by  noise,  and  in  the 
limit  of  infinitely  long  measurement  time,  will  equal  the  true  values  of  the  surface 
impedances. 

In  practice,  of  course,  measurement  times  are  not  infinitely  long.  In  this  case,  both  the 
standard  and  the  remote  reference  techniques  will  produce  estimates  of  the  surface 
impedances  that  are  different  from  the  true  values.  As  we  saw  above,  the  standard  method 
results  will  tend  to  be  smaller  in  magnitude  than  the  true  values,  while  the  remote 
reference  method  results  will  vary  both  above  and  below  the  true  value.  The  amount  of 
variation  from  the  true  value  at  any  given  frequency  depends  on  the  intensity  of  the 
electromagnetic  signal  at  that  frequency  compared  to  the  noise.  As  an  example,  Figure  1 
shows  the  magnitude  of  the  impedance  computed  from  a  short  data  record.  The  top  curve 
is  the  remote  reference  result;  the  bottom,  the  standard  result. 


Frequency  (Hz) 

Figure  1.  Standard  and  Remote  Reference  Impedance  Estimates 
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It  is  important  to  understand  that,  for  the  purposes  of  finding  the  impedance,  the  “signal” 
includes  all  electromagnetic  radiation  from  distant  sources,  whether  they  are  controlled 
sources  or  natural  random  sources  like  lightning.  Any  signal  from  a  distant  source, 
random  or  controlled,  satisfies  the  assumptions  necessary  for  the  impedance  to  be 
meaningful.  “Noise”  means  anything  that  is  not  a  propagating  EM  wave  from  a  distant 
source,  and  so  does  not  satisfy  the  assumptions  needed  for  the  impedance  to  be 
meaningful,  like  the  internal  noise  of  the  sensors,  local  electromagnetic  interference,  or 
noise  due  to  vibration  of  magnetic  sensors. 

In  Figure  1,  the  remote  reference  and  standard  method  impedances  are  different  from 
each  other  in  the  region  between  2  kHz  and  5  kHz  where  little  natural  background  signal 
is  present;  above  and  below  this  band,  the  two  impedances  are  closer  together  because  of 
the  presence  of  signals  caused  by  distant  thunderstorm  activity.  The  large  difference 
between  the  two  impedance  estimates  above  in  the  30  -  50  kHz  region  is  due  to  the 
effects  of  internal  noise  in  the  magnetic  field  sensors. 

2.1.2.  Power  Spectral  Estimation 

One  of  the  reasons  for  formulating  the  impedance  equations  as  we  have  above,  so  that  the 
entries  computed  from  the  measured  data  are  auto-  or  crosspower  spectral  densities,  is  to 
take  advantage  of  the  large  body  of  previous  work  on  accurate  estimation  of  spectral 
densities.  The  two  major  categories  of  power  spectral  estimation  techniques  are  the 
classical  methods  based  on  the  Fourier  transform,  and  the  parametric  methods,  which 
assume  that  the  signal  can  be  modeled  as  a  particular  type  of  random  process,  and 
estimate  the  model  coefficients  in  order  to  determine  the  spectral  estimate.  The  classical 
methods  are  the  oldest  and  most  widely  used,  and  although  parametric  techniques  have 
demonstrated  advantages  in  other  applications,  they  have  not  been  widely  applied  in 
processing  of  magnetotelluric  data.  At  present,  classical  methods  are  used  for  the 
processing  of  all  field  data,  but  an  examination  of  the  possible  benefits  of  parametric 
methods  is  planned  during  the  next  year. 

Because  classical  spectral  estimation  techniques  have  been  in  wide  use  for  more  than 
thirty  years,  their  proper  use  is  very  well  understood;  in  fact,  a  good  classical  spectral 
estimate  can  be  computed  in  almost  the  manner  implied  by  the  equations  given  in  the 
previous  section.  As  we  can  see,  each  entry  in  the  final  matrix  equation  is  the  product  of 
the  Fourier  transform  of  a  signal  with  the  conjugate  of  the  Fourier  transform  of  either  the 
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same  signal  (an  autopower  spectral  density),  or  a  different  signal  (a  crosspower  spectral 
density).  If  we  were  to  simply  take  the  Fourier  transform  of  the  entire  input  data  record, 
compute  its  conjugate,  and  multiply  them  together,  however,  we  would  have  a  poor 
spectral  estimate  for  two  important  reasons:  leakage  and  variability. 

Leakage  refers  to  the  spreading  of  the  energy  from  a  sharp  spectral  feature,  such  as  a 
sinusoidal  signal,  into  adjacent  frequencies,  where  it  may  obscure  other  details  of  the 
spectrum.  This  is  caused  by  the  fact  that  a  finite  data  record  is  in  effect  the  product  of 
multiplying  an  infinite  data  record  with  a  rectangular  “window”  function  that  is  one 
during  the  observation  time,  and  zero  elsewhere.  In  the  Fourier  domain,  multiplication  of 
the  data  by  a  window  is  equivalent  to  the  convolution  of  the  Fourier  transform  of  the  data 
with  the  Fourier  transform  of  the  window,  and  so  if  the  window  has  a  Fourier  transform 
that  is  spread  out  across  a  range  of  frequencies,  any  sharp  spectral  features  in  the  data  will 
also  be  spread  out.  A  rectangular  window  is  just  the  result  of  doing  nothing  to  the  input 
data,  but  there  are  many  other  windows  with  Fourier  transforms  that  are  more  compact  in 
the  frequency  domain  and  cause  much  less  spreading  in  the  computed  spectra.  All  that  is 
necessary  to  dramatically  reduce  leakage  is  to  choose  an  appropriate  window  function  and 
multiply  the  data  record  by  that  window  before  computing  the  Fourier  transform. 

To  illustrate  the  importance  of  windowing  to  reduce  leakage,  Figure  2  below  compares  an 
estimate  computed  using  a  rectangular  window  with  one  computed  using  a  low-leakage 
Hann  window,  on  an  actual  data  record  taken  with  a  local  source  operating  at  5  kHz.  The 
figure  shows  the  field  spectral  density,  which  is  the  square  root  of  the  power  spectral 
density,  as  is  more  common  in  geophysics.  Notice  that  when  the  rectangular  window  is 
used,  the  strong  5  kHz  signal  is  spread  over  a  very  large  range,  completely  obscuring 
other  spectral  features.  The  Hann  window  has  dramatically  lower  leakage,  providing  an 
accurate  estimate  of  spectral  features  near  the  powerful  5  kHz  signal. 
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Figure  2.  Proper  Choice  of  Window  is  Essential  when 
Sharp  Spectral  Features  are  Present 


The  second  problem  that  must  be  addressed  in  classical  spectral  estimation  is  the 
variability  of  the  estimate.  If  we  compute  the  windowed  Fourier  transform  of  the  entire 
input  data  record,  we  will  have  a  spectral  estimate  with  very  high  resolution,  but  also  very 
high  variability.  The  standard  deviation  of  a  spectral  estimate  computed  in  this  way  is 
approximately  equal  to  the  true  spectral  density,  so  the  variability  of  the  estimate  is  very 
large.  To  reduce  this  variability,  we  can  divide  the  input  data  record  into  shorter 
segments,  window,  and  compute  spectral  estimates  from  each  segment,  then  average  the 
segment  estimates  together  to  form  the  overall  spectral  estimate.  This  will  reduce  the 
variability  at  the  cost  of  decreasing  resolution.  In  fact,  we  can  trade  off  variability  and 
resolution  by  increasing  or  decreasing  the  number  of  segments. 

The  combination  of  windowing  and  segment  averaging  produces  a  classical  spectral 
estimation  method  with  excellent  performance.  In  fact,  the  procedure  we  have  described 
above  is  a  very  widely  used  method  first  developed  by  Welch.  Welch  also  found  that 
using  segments  that  overlap  can  further  reduce  the  variability  of  the  spectral  estimate;  in 
computation  of  estimates  on  our  experimental  data,  we  have  employed  a  Hann  window  to 
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minimize  leakage,  and  overlapped  the  data  segments  by  50%  to  minimize  the  variability 
of  the  estimate. 

As  we  have  seen,  the  major  tradeoff  in  the  computation  of  power  spectral  estimates  using 
classical  techniques  is  between  the  resolution  of  the  spectral  estimate  and  its  accuracy. 
From  a  given  amount  of  input  data,  it  is  possible  to  compute  spectral  estimates  with 
higher  resolution  and  higher  variability,  or  lower  resolution  and  lower  variability.  To 
illustrate  this  tradeoff,  we  have  computed  power  spectral  estimates  of  varying  resolutions 
from  the  same  ten-second  data  record  taken  with  a  local  source  operating  at  5  kHz. 


Figure  3.  Spectral  Estimate  for  4096  Point  FFT  Length 


The  first  estimate,  shown  in  Figure  3,  uses  a  FFT  length  of  4096  points,  resulting  in  an 
approximate  frequency  resolution  (spectral  bin  width)  of  24.4  Hz.  Although  most  of  the 
spectral  features  above  1  kHz  are  visible  in  the  spectrum,  the  resolution  is  too  low  to 
show  details  below  a  few  hundred  Hertz.  In  addition,  the  strong  peak  at  5  kHz  still  shows 
a  small  amount  of  leakage  at  its  base,  despite  the  use  of  a  low  leakage  window.  However, 
the  continuous  portions  of  the  spectrum,  where  no  sine  waves  are  present,  show  very  little 
variability;  for  example,  the  shoulder  from  6  kHz  to  10  kHz  is  very  smooth.  This 
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spectrum  is  the  average  of  slightly  more  than  500  overlapping  spectra,  reducing  the 
variability  of  the  estimated  field  spectrum  to  a  very  low  level. 


Figure  4.  Spectral  Estimate  for  16384  Point  EFT  Length 


The  second  estimate  (Figure  4)  uses  a  FFT  length  of  16384,  for  a  frequency  resolution  of 
6.1  Hz.  At  this  higher  resolution,  spectral  detail  is  now  visible  down  to  a  few  hundred 
Hertz,  and  all  of  the  sinusoidal  peaks  in  the  spectrum  are  much  better  defined.  Note  that 
as  the  resolution  increases,  the  leakage  tails  around  strong  signals  like  the  5  kHz 
transmitter  signal  become  narrower,  and  the  peaks  themselves  increase  in  height.  This 
allows  the  detection  of  smaller  peaks  like  the  cluster  visible  at  3  kHz.  However,  note  that 
the  variability  of  the  spectrum  has  also  increased;  the  shoulder  region  at  6  -  10  kHz  now 
shows  noticeable  roughness,  even  though  the  true  spectrum  is  smooth  in  this  region.  This 
estimate  is  the  average  of  approximately  130  overlapped  spectra.  Each  increase  by  a 
factor  of  four  in  the  FFT  length  reduces  the  number  of  spectra  averaged  by  a  factor  of 
four,  and  increases  the  variance  of  the  field  spectral  density  estimates  shown  by  a  factor 
of  2. 
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Figure  5.  Spectral  Estimate  for  65536  Point  FFT  Length 


The  final  example  (Figure  5)  uses  a  FFT  length  of  65536  points,  for  a  frequency 
resolution  of  1.5  Hz.  Spectral  features  are  now  visible  below  100  Hz,  and  the  sinusoidal 
peaks  have  continued  to  increase  in  height.  However,  this  spectrum  is  the  average  of  only 
31  overlapped  spectra,  and  the  variability  is  becoming  quite  large. 

We  have  seen  that  classical  power  spectral  estimation  using  a  fixed  amount  of  input  data 
involves  a  tradeoff  between  resolution  and  variability.  In  processing  the  field  data  for  this 
program,  data  records  from  10  seconds  to  100  seconds  in  length  have  been  employed. 
Typical  FFT  lengths  have  been  from  4096  to  16384  points,  resulting  in  power  spectral 
estimates  similar  to  those  in  Figures  3  and  4.  Although  these  parameters  have  yielded 
good  results,  the  interaction  between  spectral  estimate  resolution  and  variability  and  the 
performance  of  the  imaging  algorithm  has  not  yet  been  fully  explored.  A  detailed  study  of 
these  interactions  will  be  conducted  during  the  next  year. 

2.1.3.  Implementation 

The  spectral  estimation  algorithms  described  above  have  been  implemented  in  a 
processing  program  for  analysis  of  field  data.  The  program,  written  in  C,  performs  power 
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spectral  estimation  using  the  Welch  method,  applies  instrument  calibrations,  and 
produces  auto-  and  crosspower  spectral  estimates.  In  addition,  the  program  computes  the 
standard  error  (normalized  standard  deviation)  of  each  of  the  autopower  spectral 
estimates,  to  provide  an  objective  measure  of  the  quality  of  each  spectral  estimate.  All 
computations  are  made  in  32-bit  floating  point  to  ensure  accuracy,  and  using  optimized 
assembly  implementations  of  the  fast  Fourier  transform  specialized  for  real  input  data  to 
maximize  speed. 

2.2.  Imaging  Algorithms 

After  the  spectral  estimates  have  been  computed,  and  a  surface  impedance  calculated 
from  these  spectral  estimates,  the  data  are  used  to  produce  an  image  showing  the 
estimated  subsurface  conductivity  as  a  function  of  depth.  In  contrast  to  the  spectral 
estimation  problem,  which  is  very  stable,  the  inverse  imaging  problem  is  unstable,  and  a 
variety  of  methods  have  been  employed  in  the  past  to  achieve  stability  while  still 
maximizing  resolution. 

In  the  following  section,  a  new  approach  developed  by  APTI,  called  spectral 
regularization,  is  described  and  compared  with  conventional  magnetotelluric  inversion 
algorithms.  Although  the  algorithm  was  not  developed  under  this  contract,  a  basic 
understanding  of  its  operation  and  its  differences  from  conventional  magnetotelluric 
inversion  algorithms  is  important  for  understanding  the  analysis  of  experimental  data 
given  in  Section  4. 

2.2.1.  Fundamentals  of  Inverse  Imaging 

One  often  has  a  set  of  observed  data  e  arrived  at  as  the  result  of  an  experiment,  that  one  is 
attempting  to  fit  with  some  function/.  This  fit  is  defined  so  as  to  minimize  an  associated 
cost  C\f\ ,  while  also  ensuring  that  a  distortion  function  D[eJ],  measuring  the  difference 
between  the  observed  data  and  the  function  /  does  not  exceed  some  maximum  value.  In 
mathematical  notation  we  have 

minC[/]  subject  to  D[e,f]<D,  D>  0  (I) 
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where  /  e  F  defines  the  class  of  admissible  solutions.  It  is  important  to  appreciate  that 
problem  (I)  has  a  dual  that  can  be  written  as 

min D[e,f]  subject  to  C[f]<C,  C  <°°,  (II) 

feF 

in  other  words,  minimize  the  distortion  while  keeping  the  cost  from  exceeding  a  given 
value.  In  both  cases,  critical  choices  must  be  made  in  the  specification  of  F,  C\f\,  and 
D[ef\.  Some  choices  result  in  a  problem  that  has  no  solution;  in  other  cases  the  solution 
may  not  be  unique,  stable,  or  computationally  feasible  to  determine. 

Analytically,  the  use  of  Lagrange  multipliers  provides  us  with  some  basic  insight  into  the 
connections  between  (I)  and  (II).  For  any  X  >  0 ,  the  optimum  solution,  f*x ,  for  the 

unconstrained  problem  defined  by  the  Lagrange  objective 

JJf\e]  =  D[e,f]  +  XC[f], 

written  in  the  usual  way  as 

f*  =  argmin/^t/je], 

feF 

is  also  optimal  for  the  constrained  problem 

fl  =  arg  min  C[f  ]  subject  to  D[e,f]<  D[e,  f*x  ] . 

/eF 

Hence,  X >  0  indexes  a  family  of  achievable  solutions  with  the  property  that  if  D[e,fx] 
equals  D ,  then  f*x  is  an  optimal  solution  to  (I).  Emphasizing  duality,  fx  is  also  optimal 
for  the  constrained  problem  defined  according  to 

fx  =  arg  nun  D[e,  f  ]  subject  to  C[f]<  C[fx  ] , 
and  if  C[fx]  equals  C ,  then  fl  is  an  optimum  for  (II). 


15 


In  general,  there  exists  a  set  of  pairs,  {c[fl  ],  D[e,  fl  ]) ,  indexed  by  A ,  that  minimize  the 

unconstrained  Lagrangian  and  occupy  the  space  C  x D,  upon  which  all  the  extremals  are 
defined. 

2.2.2.  Conventional  Magnetotelluric  Methods  and  New  Approaches 

The  recent  literature  on  MT  inversion  focuses  on  a  type  (I)  constrained  optimization, 
defining  D[e,f]  as  a  squared  error  and  modeling  this  error  as  Gaussian  so  that  %2 
random  variables  result.  The  constraint  is  expressed  as  a  hard  constraint  of  the  form 

D[e,f]  =  D, 

and  the  Lagrangian  formulation  is  used  sequentially  to  estimate  a  limiting  operating  point, 
A*,  that  satisfies  this  constraint  while  delivering  a  minimum  cost  C[//] .  In  these 

formulations,  the  fit  is  defined  as  a  mapping  from  spatial  coordinates  to  conductivity 
values  so  that  D[e,f£]  is  of  the  form  where  o^>  0  is  the  desired 

conductivity  estimate.  This  reflects  the  fact  that  the  computation  of  the  distortion  entails 
the  forward  (numerical)  solution  of  Maxwell’s  equations.  In  addition,  these  inversion 
methods  choose  C[f]  as  a  measure  of  variability  or  roughness,  which  is  motivated  in 
general  by  Tikhonov  theory  as  a  way  to  deal  with  the  instability  of  the  inverse  problem. 

With  reference  to  the  preceding  discussion,  our  current  approach  differs  from  most  MT 
methods  in  two  key  ways.  First,  we  use  a  type  (II)  formulation  and  place  the  constraint  on 
the  cost  C[/]  while  the  distortion  D[e,f]  is  minimized.  In  this  way,  each  potential  fit  is 
subjected  to  the  same  total  cost,  while  the  optimal  fit  is  the  one  that  minimizes  the  total 
distortion.  Of  course,  the  mathematics  are  dual  and  one  can  pose  the  optimization  in 
either  way,  but  the  physical  aspects  of  the  problem  suggest  that  option  (II)  has  some 
advantages.  Second,  we  divide  the  problem  into  two  nested  optimizations,  in  which  the 
first  obtains  an  estimate  of  the  resistivity  at  each  surface  point  as  a  function  of  frequency, 
using  the  measured  field  data,  and  the  second  addresses  the  spatial  inversion  aspect.  So, 
in  fact,  we  employ  two  type  (II)  optimizations  in  order  to  obtain  the  final  result.  Having 
made  this  separation,  we  may  then  define  a  cost  for  the  first  optimization  that  follows 
directly  from  physical  considerations.  This  is  in  contrast  to  the  costs  used  in  many  MT 
algorithms  that  are  largely  chosen  without  regard  to  the  physics  of  the  problem. 
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2.2.3.  Spectral  Regularization 


To  be  specific,  the  method  employed  here  differs  from  earlier  methods  in  that 
fundamental  physical  constraints,  namely  the  Kramers-Kronig  relations  that  arise  in 
dispersive  wave  theory,  dictate  the  form  of  the  C[f]  used  in  the  first  optimization.  This 
approach  is  unique  in  that  it  employs  a  filtering  of  the  raw  data,  referred  to  as  spectral 
regularization ,  before  any  inversion  is  performed.  This  type  of  pre-inversion  filtering 
yields  estimates  that  can  then  be  either  jointly  processed  or  post-processed  using 
conventional  MT  inversion  techniques  to  achieve  resistivity  imaging  with  higher 
resolution.  Our  approach  also  differs  in  two  other  key  aspects,  but  these  aspects  are  better 
discussed  after  the  basic  principles  have  been  explained. 

The  general  form  of  the  optimization  employed  is  most  easily  expressed  in  a  two-step 
format.  The  first  step  of  the  optimization  involves  minimizing,  for  a  given  coherence 
threshold  k,  and  spectral  regularization  parameter  Av,  a  vertical  sounding  objective 

JVI X.^vKPc)  with  respect  to  pc ,  where  Jv  is  given  as 


1) 


[K,^\pc](pc)  =  \CK(f) 


.  Pcif) 

logm 


4f  +  ^v| 


d  log  Pe  CjO 
dlogf  ' 


Here,  the  input  pc  is  the  complex  resistivity  defined  as 


Pdf’*) 


Z\f\x) 


where  Z(f;x)  is  an  estimate  of  the  surface  impedance  Z(f  ;x)  computed  from  the 
measured  data  at  each  frequency  /  and  position  x.  Finally,  CK(f;x )  is  the  thresholded 
coherence  matrix  defined  as 

CK(fix)  =  l{C(f;x)>K}  0<k<1, 


where  K  is  the  threshold  and  C{f\x)  is  the  estimated  coherence.  The  coherence  at  any 
given  frequency  between  the  measured  electric  and  magnetic  field  associated  with  a 
particular  mode  is  a  measure  of  how  consistent  the  actual  measured  data  are  with  the 
underlying  assumption  of  plane  wave  excitation.  Accepting  only  data  with  a  coherence 
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above  a  certain  threshold  is  a  simple  but  effective  means  of  preventing  input  data  severely 
corrupted  by  noise  from  affecting  the  final  solution. 

If  we  define  the  solution  to  the  first  (vertical)  optimization  as 

Pi  =  argmin  Jv  [K,Xv]{pc) , 

Pc 


then  we  may  define  the  second  (horizontal)  optimization  objective  as 


2) 


Jbfc»|p;](pe)=j  pcw-p;w!*+aJ 


4>c(*) 

fdx 


dx . 


We  have  employed  the  shorthand  notations  pc(f)  and  pvc(x)  for  pc(f;x )  and  pl(x;f) , 
respectively,  to  keep  the  notation  simple.  If  we  define  the  nested  objective 


J[K,Av,Ah]  =  Jl 


A 


argmin  J 

Pc 


[K,A  PcK Pc) 


then  the  overall  optimal  solution  may  be  concisely  expressed  as 

Pc  =  arg  min  J  [k,A  ’A  KPc )  • 

Pc 

We  first  note  that  the  Jv  (vertical)  spectral  regularization  objective  is  central  to  the  high 
performance  of  this  technique.  This  first  optimization  process  may  either  be  considered  as 
data  filtering  in  preparation  for  the  inversion  step,  or  included  in  the  inversion  objective 
itself. 

Here,  we  may  point  out  the  two  additional  aspects  of  this  approach  mentioned  above  that 
appear  to  also  differ  from  current  MT  algorithms.  Referring  to  Jv,  we  see  that  this 

objective,  both  its  familiar  Hilbert  space  Z?  error  integral  and  its  Sobolev  smoothing 
integral,  are  defined  in  terms  of  the  logarithm  of  the  complex  resistivity.  To  simplify 
further  discussion  and  notation,  we  will  define  A  =  log  pc  which  allows  us  to  write  the 

spectral  regularization  objective  in  the  compact  form 
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Jv[«r, Xi A](A)  =  JcJa- a| 2df  +  xj  I /  Af  df  . 


Secondly,  we  note  that  we  have  continued  to  employ  the  coherence  mask  CK  in  the 
objective  definition  even  in  cases  where  the  input  data,  A ,  have  been  derived  using  the 
remote  reference  technique.  We  next  must  turn  the  discussion  to  how  the  spectral 
regularization  objective  Jv  derives  its  form  from  basic  physical  principles. 

2.2.4.  A  Physical  Interpretation  of  Spectral  Regularization 

Consider  the  estimation  of  a  one-dimensional  ground  conductivity  o (z) 1  over  some  range 
of  depth,  z,  given  the  impedance  Z0(/)  =  E0(f)/H0(f )  measured  on  the  ground  surface 
over  some  range  of  frequency.  The  mapping  Mf  that  takes  of)  into  Z0(f)  is  strongly 
smoothing.  For  example,  Z0(-)  is  a  smooth  function  of  frequency  even  when  cr(  )  consists 
of  discrete  layers  of  markedly  different  conductivities.  Hence,  Mf  can  be  described  as  a 

low  pass  operator,  since  rapid  changes  in  &0) ,  those  with  substantial  high  spatial 
frequency  content,  are  smoothed  out.  This  suggests  that  the  inverse  operator  M~x  z>  0, 
must  have  the  kind  of  high  pass  character  to  produce  a  &{■)  that  changes  significantly 
with  depth  from  a  relatively  smooth  Z0(-).  Practically  speaking,  M”1  is  not  a  stable 
inverse  in  the  sense  that  small  changes  in  Z0(-)  must  invert  to  large  changes  in  cr(-). 
Unfortunately,  when  Z0(  )  is  computed  from  measured  data  for  E0  and  H0 ,  the  measured 
value  of  Z0(-)  has  errors  due  to  noise  in  the  measured  fields.  As  a  result,  any  noise 

reduction  stabilization  (NRS)  entails  the  problem  of  minimizing  errors  due  to  noise 
without  removing  small  changes  in  Z0(-)  due  to  fine  details  in  a{-). 

NRS  may  be  directly  applied  to  the  measured  data.  Commonly  used  techniques  in 
magnetotellurics  involving  NRS  on  the  input  data  employ  coherence  relationships  and 
remote  reference  processing,  the  use  of  spectral  whitening,  and  j2  filtering  constraints. 
Other  techniques  are  also  used  to  accomplish  NRS  at  the  output,  for  example,  by  using 
any  known  geologic  information  to  help  constrain  the  solution.  More  general  algorithms 
of  this  type  penalize  rapid  variations  in  the  estimate  of  <x(-)  and  require  increasing 
smoothness  with  increasing  depth.  These  algorithms  are  based  on  the  difficult  to  quantify 
but  qualitatively  obvious  truth  that  resolution  degrades  with  depth.  The  difficulty  arises 

1  We  adopt  the  convention  that  <J  will  refer  to  conductivities  that  vary  only  as  a  function  of  space, 
a  -  a(r) ,  while  1 ] p  will  refer  only  to  conductivities  that  vary  with  frequency,  i.e.,  apparent  resistivities 
of  the  form,  p=  p(f). 
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not  only  because  Mz  1  is  high  pass,  but  also  because  it  is  nonlinear,  and  furthermore 
nonlocal,  involving  all  of  Z0(-)  in  principle  for  each  z>  0. 

The  NRS  approach  represented  by  the  spectral  regularization  objective  Jv  penalizes  lack 
of  smoothness  due  to  noise  at  the  input.  Put  simply,  we  employ  a  smoothness  constraint 
on  an  input  observable  which,  in  the  noise-free  case,  is  known  to  be  very  smooth.  This  is 
in  contrast  to  applying  a  direct  spatial  smoothness  constraint  on  cr(-)  which  may  be 
significantly  non-smooth.  Indeed,  when  attempting  to  detect  manmade  structures,  the 
regions  where  <r(-)  is  not  smooth  are  of  paramount  interest.  In  addition,  the  input 
smoothness  constraint  we  describe  is  derived  from  the  physics  of  the  problem,  whereas 
output  smoothness  constraints  are  typically  ad  hoc,  given  the  difficulties  in  the  direct 
analysis  of  the  nonlinear  functional  operator  M~\  Nevertheless,  one  must  still  eventually 

deal  with  output  constraints  in  order  to  address  the  destabilization  implicit  in  the  forward 
smoothing  posed  by  Mf  when  the  data  are  noisy.  The  method  we  discuss  may  be  viewed 

as  taking  some  of  the  smoothing  responsibility  off  this  last  step,  and  therefore  helping  to 
maintain  maximum  resolution. 

2.2.5.  Complex  Apparent  Resistivity 

A  further  difference  appears  in  the  choice  of  the  quantity  computed  by  the  regularization 
algorithm.  The  well  known  apparent  resistivity  has  the  form 


This  contains  no  phase  information  since  the  absolute  value  of  the  impedance  has  been 
taken.  Often,  in  standard  magnetotelluric  analysis,  the  phase  of  Z  is  used  in  conjunction 
with  pa  so  that  this  phase  information  is  not  lost.  This  approach  leads  to  awkward 
calculations  since  the  pa  and  phase  functions  must  be  manipulated  separately.  A  more 
practical  definition,  which  allows  the  field  properties  of  the  complex  numbers  to  take  care 
of  both  the  magnitude  and  phase  simultaneously  during  calculation,  is  to  define 


which  we  call  the  complex  resistivity.  In  the  noise-free  case,  the  electric  field  leads  the 
magnetic  field  by  45°  in  a  resistive  half-space  so  that  argZ2  =n\ 2  and  therefore#, 

although  defined  in  terms  of  complex  quantities,  delivers  a  real  number  as  it  should  in 
this  simple  situation. 

Use  of  this  complex  resistivity  makes  it  straightforward  to  derive  the  basic  physical 
relationship  that  is  exploited  in  spectral  regularization.  For  the  one-dimensional  problem, 
the  magnitude  and  phase  of  the  complex  resistivity  fie  are  related  by  the  well  known 
Kramers-Kronig  dispersion  relations.  For  example,  the  phase  can  be  expressed  in  terms  of 
the  logarithm  of  the  magnitude  of  the  complex  resistivity  as 

4>c(/)  =  fjri0^(<p)|y2^2  > 

where  <7  is  taken  as  the  "background"  conductivity  for  which  <7  =  <7(°°)  .  This  relation 
holds  without  reference  to  the  good  conductor  approximation  and  so  is  valid  at  any 
frequency.  In  the  case  of  a  constant  half-space,  Pc{f)  =  V<7  for  all  / ,  not  merely  in  the 
limit  for  well  defined  cfz),  and  so  log|cpc|  =  0,  implying  Zfidf)  =  0,  i.e.,  pdf)  is 
real,  as  mentioned  above. 


A  logarithmic  change  of  variables  converts  the  above  expression  into  a  convolution  with 
an  odd  kernel.  Integration  by  parts  therefore  results  in  an  even  kernel,  one  having  very 
rapid  decay,  which  when  approximated  using  a  Dirac  ^-function  results  in 


7tdlog|pc(/)| 
4  dlogf 


This  approximate  expression  is  well  known  in  the  magnetotelluric  community,  and 
provides  a  good  approximation  for  cases  with  one-dimensional  conductivity  variation.  By 
using  this  expression,  we  can  write 


dA 


d  log  f  dlogf 


[Log|pc|  +  ;Zpc] 


_  2  .  dZPc 

n  c  ^  dlogf 


with  the  result 
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In  other  words,  minimizing  the  prior  component  of  the  Jv  objective,  as  suggested  by  the 
Kramers-Kronig  relation,  may  be  viewed  as  minimizing  a  weighted  sum  of  the  squares  of 
the  differences  between  the  observed  phase  and  group  phase  and  those  that  would  be 
observed  for  a  space  with  only  a  one-dimensional  conductivity  variation. 

This  approximation  can  also  be  expressed  in  an  equivalent  form 


/A'  =  / 


Pc 


which  is  particularly  convenient  when  the  complex  resistivity  parameterization  is  used. 
Like  conventional  regularizers  based  on  derivatives  of  the  apparent  resistivity,  this 
method  stabilizes  the  initial  ill-posed  problem.  However,  unlike  conventional 
regularizers,  its  form  is  based  on  a  reasonable  physical  approximation  rather  than  an  ad 
hoc  choice  made  for  computational  convenience. 

2.2.6.  Summary  of  Imaging  Method 

The  expressions  derived  above  are  of  course  approximate,  but  they  can  be  used  to  address 
the  input  NRS  problem  in  a  systematic  fashion  based  on  physical  principles.  The 
convolution  expressions  (not  shown  here)  can  themselves  be  employed  as  constraints 
without  resort  to  approximation,  other  than  in  a  modeling  sense,  but  this  approach  has 
some  difficulties  that  have  yet  to  be  addressed.  We  have  described  the  reasoning  leading 
to  an  optimization  objective  of  the  form 

J v  [K,X](A)  =  J  cK  |a  -  a| 2df  +  xj  1/  Af  df  . 


Minimization  of  the  above  form  results  in  an  optimal  A  that  is  conditioned  on  an 
approximate  half-space  model.  Hence,  solution  anomalies  may  then  be  interpreted  as  due 
to  physical  departures  from  this  half-space  model.  Several  generalizations  of  the  basic 
approach  follow  immediately. 


The  basic  principle  is  clear:  before  any  inversion  proper,  it  is  useful  to  filter  the  data  so 
that  it  more  closely  respects  the  Kramers-Kronig  relation.  This  physically  justifiable 
smoothing  process  is  applied  to  noisy  data  that  is  known  to  be  smooth  when  no  noise  is 
present.  The  form  of  the  filter  is  simple:  more  smoothing  is  applied  at  the  higher 
frequencies  than  the  lower  ones.  The  above  Lagrangian  form  can  of  course  be  posed  using 
different  metrics.  From  the  computational  point  of  view  each  of  the  above  optimizations 
can  be  implemented  using  fast  and  stable  signal  processing  techniques. 
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3.  Electromagnetic  Modeling  and  Inversion 


Previous  theoretical  and  computational  investigations  of  underground  conductivity 
imaging  using  low  frequency  methods  have  been  performed  mostly  for  geological 
structures;  there  has  been  very  little  work  performed  in  modeling  underground  targets  of 
military  interest.  Consequently,  this  contract  includes  a  significant  modeling  effort 
specifically  aimed  at  determining  the  characteristics  of  man-made  targets.  The  results  of 
this  effort  are  being  used  to  plan  experiments,  evaluate  and  improve  the  performance  of 
signal  processing  and  inversion  algorithms,  choose  source  signals,  and  predict 
observables  associated  with  complex  underground  structures. 

To  perform  this  work,  APTI  has  employed  a  set  of  standard  geophysical  simulation  and 
inversion  programs,  with  parameter  settings  properly  adapted  to  the  characteristics  of 
man-made  targets,  to  study  the  response  of  underground  structures  of  military 
significance.  The  results  of  these  efforts  to  date  are  described  in  the  following  sections. 

In  Section  3.1,  we  present  the  results  of  two-dimensional  forward  simulations  of  an 
underground  tunnel,  both  in  a  half-space  with  uniform  conductivity  and  in  the  presence  of 
thin  conductive  surface  layers.  These  forward  simulations  were  performed  using  the  finite 
element  code  PW2D,  developed  by  Wannamaker  and  Stodt1.  Several  cases  with  varying 
ratios  of  tunnel  depth  to  tunnel  size  were  investigated.  In  Section  3.2,  we  present  the 
results  of  inverting  this  simulated  data  with  a  standard  geophysical  inversion  code,  the 
Rapid  Relaxation  Inverse  program  developed  by  Smith  and  Booker2,  and  an  analysis  of 
the  results. 


1  P.E.  Wannamaker  and  J.A.  Stodt,  A  Stable  Finite  Element  Solution  for  Two-Dimensional  Magnetotelluric 
Modeling.  Geophys.  J.R.  Astr.  Soc.  88  ,  pp.  277-296  (1987). 

2  J-T.  Smith  and  J.R.  Booker,  Rapid  Inversion  of  Two-  and  Three-Dimensional  Magnetotelluric  Data.  J. 
Geophys.  Res.  96  ,  No.  B3,  pp.  3905-3922  (1991). 
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3.1.  Forward  Modeling 


An  example  of  the  models  used  in  the  forward  simulations  is  shown  in  Figure  6.  This 
model  contains  a  uniform  half-space  with  resistivity  100  Ohm-m,  and  a  15  xl5  m  tunnel 
at  a  depth  of  50  m.  From  this  base  situation  variations  in  tunnel  size,  as  well  as  the  effects 
of  conducting  surface  layers,  were  examined.  The  resistivity  of  the  tunnel  interior  was 
taken  as  104  Ohm-m.  In  all  cases  discussed  here  the  TM  mode  (surface  E  field  in  the  .re¬ 
direction)  was  used;  this  is  the  most  significant  mode  for  the  resistive  tunnels  considered 
here.  The  frequencies  used  in  the  simulations  spanned  the  range  from  100  Hz  to  100  kHz. 
For  forward  model  creation,  the  UNIX  based  graphical  user  interface  MTWorks  was 
used.  Further  calculations  were  performed  by  exporting  the  results  from  MTWorks  and 
importing  them  into  Matlab.  Computations  were  performed  on  a  Sun  Ultra  140  with  320 
Megabytes  of  RAM;  computation  time  for  a  typical  forward  model  with  10,000  elements 
was  of  the  order  of  two  to  three  hours  for  45  frequencies. 

3.1.1.  Uniform  Half  Space 

The  objective  of  the  first  set  of  simulations,  whose  results  are  shown  in  Figures  7-10,  was 
to  investigate  the  effect  of  the  depth-to-size  ratio  on  the  field  perturbation  at  the  surface, 
as  a  function  of  frequency  and  horizontal  distance  from  the  tunnel  center.  The  response  to 
square  tunnels  with  sizes  5  meters  (Figure  7),  10  meters  (Figure  8),  15  meters  (Figure  9) 
and  25  meters  (Figure  10)  was  studied.  These  correspond  to  depth-to-size  ratios  of  2,  3.3, 
5,  and  10.  The  right-hand  sides  of  Figures  7-10  show  the  apparent  resistivities  while  the 
left-hand  sides  show  the  phase  of  the  impedances  as  a  function  of  the  frequency  and 
horizontal  position.  As  a  benchmark,  note  that  the  skin  depth  at  8  kHz  is  approximately 
56  m,  close  to  the  tunnel  depth;  we  would  expect  the  tunnel  response  to  be  most  evident 
near  that  frequency. 
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Figure  6.  Example  Finite  Element  Model  for  Forward  Simulation. 
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Figure  7.  Forward  Simulation  Results  for  5m  Tunnel  50m  Deep. 
Background  is  100  ohm-m. 
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4.519e*01 


Background  is  100  ohm-m. 


It  is  apparent  that  self-similar  patterns,  whose  intensity  falls  with  the  area  of  the  tunnel, 
are  present.  Concentrating  on  the  features  of  Figure  7,  we  note  that  the  maximum  phase 
lag  (approximately  4  degrees)  occurs  at  a  frequency  corresponding  to  a  true  depth  to  skin 
depth  ratio  of  approximately  0.8,  consistent  with  Bostick’s  analysis.  On  the  other  hand, 
the  maximum  of  the  apparent  resistivity  (130  Ohm-m)  appears  at  frequencies  well  below 
500  Hz  and  has  the  expected  form  of  an  electrostatic  umbra.  The  field  perturbations  due 
to  the  presence  of  the  tunnel  extend  to  lateral  distances  of  the  order  of  100  m.  Finally,  the 
sign  of  the  slope  in  both  phase  and  resistivity  changes  at  lateral  distances  of  the  order  of 
60  m.  Of  utmost  importance  to  our  inversion  procedure,  to  be  discussed  later,  is  the 
correlation  of  the  phase  and  apparent  resistivity  data,  which  under  half-space  conditions 
should  satisfy  the  Kramers-Kronig  relations.  From  Figures  7-10  we  see  that  over  the 
entire  frequency  region  below  20  kHz,  the  relationship  between  these  quantities  is 
different  from  the  constant  expected  for  half-space.  This  can  be  used  as  a  robust 
constraint  in  inversions  involving  localized  anomalies  due  to  the  presence  of  resistive  or 
conducting  tunnels. 

The  effect  of  the  tunnel  size  on  the  measurements  becomes  clearer  by  referring  to  Figures 
11  and  12,  which  show  the  differences  between  the  apparent  resistivity  and  phase  directly 
above  the  tunnel  and  at  a  horizontal  offset  of  50  meters.  Both  quantities  scale  as  the 
inverse  of  the  area  occupied  by  the  tunnel,  as  expected  on  the  basis  of  theoretical 
considerations. 

Figures  13  and  14  show  the  differences  of  measurements  taken  between  stations  located 
at  25  and  50  meter  horizontal  offsets.  These,  in  conjunction  with  the  results  shown  in 
Figures  11  and  12,  allow  us  to  characterize  the  approach  to  half-space  response  with 
distance  from  the  tunnel.  It  can  be  seen  that  the  response  is  self-similar  with  respect  to 
frequency,  but  the  amplitude  of  the  deviation  scales  with  the  inverse  square  of  the  lateral 
distance. 
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Resistivity  (Percent  Change) 


5  m  Tunnel 

dat2rho5mxyz  :  0  m  -  50  m  Station 


10  m  Tunnel 

dat2rhol0mxyz  :  0  m  -  50  m  Station 


Frequency  (Hz) 
25  m  Tunnel 


Figure  11.  Vertical  Cuts  of  Data  in  Figures  7  -  10.  Change  in 

Apparent  Resistivity  Directly  Above  Tunnel  vs.  50m  Away 
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Apparent  Resistivity  (Percent  Change) 


5  m  Tunnel 

dat2rho5mxyz  :  25  m  -  50  m  Station 


dat2rhol5mxyz  :  25  m  -  50  m  Station 


10  m  Tunnel 


dat2rhol0mxyz  :  25  m  -  50  m  Station 


25  m  Tunnel 

dat2rho25mxyz  :  25  m  -  50  m  Station 


Figure  13.  Vertical  Cuts  of  Data  in  Figures  7-10.  Change  in 
Apparent  Resistivity  from  25  -  50m  Horizontal  Distance 
from  Tunnel. 
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hange) 


5  m  Tunnel 

dat2phase5mxyz :  25  m  -  50  m  Station 


Frequency  (Hz) 

15  m  Tunnel 

dat2phasel5mxyz  :  25  m  -  50  m  Station 


10  m  Tunnel 


dat2phasel0mxyz :  25  m  -  50  m  Station 


25  m  Tunnel 

dat2phase25mxyz  :  25  m  -  50  m  Station 


Figure  14.  Vertical  Cuts  of  Data  in  Figures  7  - 10.  Change  in 
Phase  of  Impedance  from  25  -  50m  Horizontal  Distance 
from  Tunnel. 


3.1.2.  Effect  of  Conducting  Surface  Layers 


We  have  found  that  in  all  of  our  experimental  studies,  including  imaging  of  the  Silver 
Fox  and  the  Golden  Zone  mines,  the  effects  of  surface  layers  with  higher  conductivity 
were  important  in  estimating  the  tunnel  depth.  A  set  of  simulations  was  performed  to 
examine  the  effect  of  such  layers.  The  results  shown  below  were  conducted  for  a  tunnel 
and  background  similar  to  the  actual  Golden  Zone  survey  line  1  data,  whose  imaging  is 
discussed  in  Section  4  of  the  report.  The  tunnel  depth  in  the  simulations  was  28  m,  its 
size  2.5  x  3.0  m,  and  the  background  conductivity  200  Ohm-m.  Three  simulated  cases  are 
shown  here.  Figure  15  shows  the  simulation  results  in  the  absence  of  a  conducting  layer, 
while  Figures  16  and  17  show  results  in  the  presence  of  a  conducting  layer  with 
resistivity  20  Ohm-m  and  thickness  one  meter,  and  resistivity  10  Ohm-m  and  thickness 
five  meters,  respectively. 

As  expected,  the  benchmark  case  shown  in  Figure  15  is  essentially  similar  to  the  results 
of  Section  3.1.1,  once  the  results  have  been  scaled  for  the  changes  in  conductivity,  depth 
and  size.  The  maximum  phase  response  is  in  the  vicinity  of  20  kHz.  Figure  16  shows 
clearly  the  presence  of  a  poorly  resolved  top  conducting  layer,  as  expected.  However,  the 
large  phase  shifts  induced  by  the  layer  obscure  the  anomalies  caused  by  the  presence  of 
the  tunnel.  In  Figure  16,  the  tunnel  appears  as  a  small  anomaly  near  1  kHz.  The 
conducting  top  layer  acts  as  a  low-pass  filter,  reducing  the  high  frequency  response. 
Notice,  that  contrary  to  the  uniform  half-space  case,  when  a  thin  conducting  layer  is 
present  both  the  phase  and  resistivity  anomalies  have  the  same  frequency  dependence. 
Finally,  Figure  17  shows  that  a  conducting  layer  with  5  m  thickness  and  10  Ohm-m 
resistivity  completely  masks  the  tunnel,  at  least  within  the  frequency  range  used  here.  It  is 
clear  from  the  above  results  that  surface  conducting  layers  make  detection  of  tunnels 
more  difficult. 
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Figure  15.  Forward  Simulation  Results  for  a  Tunnel  2.5m  high, 
3.0m  wide,  28m  deep.  Background  is  200  ohm-m. 
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igure  16.  Conducting  Top  Layer  Problem.  Forward  Simulation  Results  for  a 
Tunnel  2.5m  high,  3.0m  wide,  28m  deep.  Background  is  200  ohm-m.  Top 
Layer  is  lm  thick  and  20  ohm-m. 
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Figure  17.  Conducting  Top  Layer  Problem.  Forward  Simulation  Results  for  a 

Tunnel  2.5m  high,  3.0m  wide,  28m  deep.  Background  is  200  ohm-m.  Top 
Layer  is  5m  thick  and  10  ohm-m. 


Apparent  Resistivity  (Percent  Change)  Apparent  Resistivity  (Percent  Change) 


No  Conducting  Top  Layer 

goIdzoneTunOOmlayOOohmRhoxyz  :  0  m  -  50  m  Station 


5  m  Thick  10  ohm-m  Conducting  Top  Layer 
goldzoneTun05mlay  lOohmRhoxyz  :  0  m  -  50  m  Station 


1  m  Thick  20  ohm-m  Conducting  Top  Layer 
goldzoneTun01mlay20ohmRhoxyz  :  0  m  -  50  m  Station 


Figure  18.  Conducting  Top  Layer  Problem.  Vertical  Cuts  of  Data  in 
Figures  15  -  17.  Effect  of  Conducting  Layers  on  Apparent 
Resistivity. 
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Phase  (Percent  Change)  Phase  (Percent  Change) 


No  Conducting  Top  Layer 

goldzoneT unOOmlayOOohmPhasexy z  :  0  m  -  50  m  Station 


5  m  Thick  10  ohm-m  Conducting  Top  Layer 
goldzoneTun05mlay!0ohmPhasexyz  :  0  m  -  50  m  Station 


1  m  Thick  20  ohm-m  Conducting  Top  Layer 
goldzoneTun01mlay20ohmPhasexyz :  0  m  -  50  m  Station 


Figure  19.  Conducting  Top  Layer  Problem.  Vertical  Cuts  of  Data  in 
Figures  15  -  17.  Effect  of  Conducting  Layers  on  Phase  of 
Impendance. 
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Figures  18  and  19  show  the  difference  in  apparent  resistivity  and  phase  between 
measurements  directly  above  the  tunnel  and  50  m  away.  It  can  be  seen  that  the  main 
effect  of  the  conducting  layer  is  to  shift  the  maximum  response  towards  the  lower 
frequencies,  while  the  amplitude  of  the  response  is  reduced  by  a  small  factor.  The 
presence  of  the  tunnel  can  still  be  detected  in  both  the  apparent  resistivity  and  the  phase 
curves.  The  results  shown  above  give  further  credence  to  inversion  techniques  using 
Kramers-Kronig  relations  as  robust  indicators  of  the  presence  of  tunnel  induced 
anomalies,  even  in  the  presence  of  thin  conducting  layers  They,  furthermore,  indicate  that 
use  of  an  equivalent  half-layer  can  be  useful  in  inversions.  This  issue  will  be  addressed  in 
Section  3.3.3. 

3.2.  Inversion  of  Modeled  Data 

To  explore  the  limits  of  conventional  techniques  used  in  magnetotelluric  exploration  for 
imaging  resistive  underground  structures  of  military  interest  (tunnels  several  diameters 
deep),  we  have  performed  inversions  of  the  noise-free  forward  simulations  presented 
above.  To  establish  a  benchmark  for  inversions  of  such  structures,  we  chose  a  widely 
used  geophysical  inversion  code,  the  RRI  inversion  code.  For  the  inversions,  we  used 
forward  simulated  data  for  underground  tunnels  in  uniform  backgrounds  with  depth-to- 
tunnel  size  ratios  between  3  and  10. 

The  results  discussed  below  show  that  with  simulated  noise-free  forward  data  accurate  to 
0.1  percent,  RRI  can  produce  a  discernible  image  of  a  tunnel  with  a  depth-to-size  ratio  of 
10.  Although  the  contrast  between  the  image  and  the  background  is  small,  and  imaging  of 
a  tunnel  with  this  depth-to-size  ratio  is  probably  not  possible  using  RRI  on  actual 
measured  (not  simulated)  data,  this  result  is  itself  unexpected.  We  have  also  found  that 
the  RRI  has  poor  resolution  on  the  small  targets  used,  even  when  they  are  large  enough  to 
be  clearly  visible  in  the  inversions. 
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In  addition,  we  have  studied  how  the  sampling  density  of  the  surface  measurements  and 
the  presence  of  conducting  surface  layers  affect  the  resolution  of  the  inversion.  We  found 
that  the  top  conducting  layer  introduces  a  second  conductivity  contrast  scale  into  the  RRI 
inversion  that  can  mask  the  presence  of  tunnels  even  though  the  tunnel  is  still  detectable 
in  the  forward  model  data  by  comparing  soundings  above  the  tunnel  with  those  far  from 
the  tunnel.  To  solve  this  multiple  scale  problem,  we  discuss  an  approach  -  mapping  two 
layer  tunnel  inversions  into  equivalent  uniform  half-space  tunnel  inversions  -  that  may 
improve  detection  and  imaging  sensitivity.  Selected  examples  of  the  RRI  inversion 
studies  are  presented  below. 

3.2.1.  Resolution  of  RRI  Inversions 

We  consider  first  the  RRI  inversion  of  a  forward  simulation  for  a  15  x  15  meter  tunnel  50 
meters  deep.  The  model  used  as  input  for  the  forward  simulation  is  shown  in  Figure  6, 
and  produces  the  forward  data  shown  in  Figure  8.  Figure  20  shows  the  result  of  an  RRI 
inversion  of  these  data.  For  this  inversion,  the  forward  data  were  sampled  every  5  m 
between  -80  m  to  +80  m.  Data  points  at  -100  m,  -90  m,  90  m,  and  100  m  were  also 
included.  The  inversion  mesh  extended  horizontally  from  -100  m  to  100  m  at  2.5  m 
intervals  and  from  zero  to  100  m  in  depth  at  5.0  m  intervals.  The  initial  model  for  the  RRI 
inversions  was  a  uniform  half-space  of  100  Ohm-m  resistivity. 

Three  tunnel  image  distortion  effects  are  apparent  in  Figure  20.  First,  the  apparent  tunnel 
depth  is  90  m  rather  than  50  m.  Second,  the  horizontal  size,  as  measured  by  the  width  of 
the  highest  resistivity  region  of  the  image,  is  25  m  —  over  60  percent  larger  than  the 
model  tunnel  width.  Third,  the  vertical  size  in  the  image  is  almost  five  times  as  large  as 
the  model  tunnel  height.  On  the  other  hand,  we  note  that  while  the  maximum  numerical 
contrast  in  the  forward  model  data  is  a  10  percent  increase  in  apparent  resistivity  above 
background  at  100  Hz,  the  inversion  produces  a  peak  increase  of  19  percent  above  the 
background  resistivity. 
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Figure  20.  RRI  Inversion  of  Forward  Data  from  Figure  9. 
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Figure  21 .  RRI  Inversion  of  Forward  Data  from  Figure  8. 
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Figure  22.  Finer  Vertical  Mesh  (Relative  to  Figure  21)  for 
RRI  Inversion  of  Forward  Data  from  Figure  8. 
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Figure  23.  RRI  Inversion  of  Forward  Data  from  Figure  7. 
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Figures  21-23  show  RRI  inversions  of  smaller  tunnels,  and  illustrate  the  effect  of  changes 
in  inversion  parameters  on  the  images  produced.  The  first  example,  Figure  21,  shows  a 
10-  meter  tunnel  inverted  using  the  same  parameters  as  Figure  20.  Notice  that  the 
resolution  is  poorer  and  the  contrast  between  the  tunnel  and  the  background  has  been 
reduced. 

To  improve  the  resolution,  the  vertical  mesh  spacing  was  decreased  by  a  factor  of  2, 
resulting  in  a  vertical  meshing  interval  of  2.5  m  from  zero  to  100  m  depth.  The  result  of 
the  inversion  is  shown  in  Figure  22.  The  increased  vertical  resolution  allows  RRI  to 
compute  more  accurate  forward  data.  The  effect  of  this  increased  accuracy  on  the  misfit 
calculation,  and  hence  resultant  inversions,  depends  on  the  size  of  the  object  of  interest 
relative  to  the  size  of  the  mesh.  For  the  15  x  15  m  tunnel,  a  vertical  mesh  of  5  m  is  one- 
third  of  a  tunnel  size;  for  the  10  x  10  m  tunnel,  a  vertical  mesh  of  5  m  is  one-half  a  tunnel 
size.  This  is  rather  coarse  for  accurate  calculations  of  the  fields  in  the  tunnel  vicinity. 
Changing  from  a  5  m  to  a  2.5  m  mesh  improves  the  accuracy  of  the  forward  computations 
within  RRI  and  results  in  a  large  difference  in  the  resulting  inversion  image. 

Figure  23  shows  the  results  of  RRI  inversion  for  the  5  x  5  meter  tunnel  at  50  meters  depth 
shown  in  Figure  10.  This  case  corresponds  to  a  depth-to-size  ratio  of  10.  The  inversion 
mesh  is  identical  to  that  used  in  the  previous  case.  However,  the  horizontal  sampling 
density  at  the  surface  in  the  region  between  -20  m  to  20  m  has  been  increased  from  every 
5.0  m  to  every  2.5  m.  The  increased  sampling  frequency  produced  an  RRI  image  superior 
to  the  previous  depth-to-ratio  case.  A  comprehensive  study  of  the  tradeoffs  between  mesh 
resolution,  horizontal  sampling  distance,  accuracy  and  computational  speed  is  in 
progress. 

The  results  in  Figure  23  show  that,  in  the  noise-free  case  and  in  a  uniform  background, 
starting  with  a  forward  simulation  of  a  tunnel  with  a  depth-to-size  ratio  of  10,  RRI  will 
produce  a  discernible  image  of  the  tunnel.  To  achieve  this,  it  was  first  necessary  to  ensure 
that  in  the  initial  forward  simulation,  numerical  errors  due  to  meshing  were  kept  small 
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(below  0.1%)  compared  to  tunnel  induced  field  perturbations  at  the  surface  of  the  earth 
model  (on  the  order  of  0.2  to  1%).  Second,  it  is  also  necessary  that  the  RRI  inversion 
mesh  be  set  fine  enough  to  do  accurate  forward  calculations  for  the  objects  of  interest. 
Third,  the  spatial  sampling  density  at  the  surface  must  be  high  enough  (at  least  every  5  m) 
to  accommodate  expected  variations  of  the  field  for  the  overall  object  size  and  depth.  The 
scale  of  this  sampling  density  should  be  a  fraction  of  the  depths  of  interest.  Fourth,  by 
making  sure  the  forward  data  was  accurate  to  0.1%,  it  is  possible  to  set  the  target  RRI 
misfit  (RMS  error)  to  0.1%.  RRI  ignores  fluctuations  less  than  the  target  misfit,  so  a  1% 
misfit  would  miss  the  majority  of  perturbations  induced  by  a  tunnel  with  a  depth  to  size 
ratio  of  10. 

From  the  viewpoint  of  geophysical  experiments  that  often  have  errors  in  field  data  on  the 
order  of  5%,  objects  that  give  rise  to  perturbations  much  smaller  than  5%  have  not 
merited  a  great  deal  of  attention.  However,  from  the  viewpoint  of  testing  the  limits  of 
detectability,  we  see  that  underground  structures  with  depth-to-size  ratios  of  10  are  within 
the  theoretical  limit  of  detection  even  with  current  algorithms  such  as  RRI. 

While  the  tunnel  in  the  forward  models  is  two  orders  of  magnitude  more  resistive  than  the 
background  (104  Ohm-m  tunnel  vs.  100  Ohm-m  background),  the  maximum  resistivity 
contrast  between  the  tunnel  image  and  tunnel  background  in  the  inversion  is  much  less. 
Further,  this  maximum  contrast  decreases  with  decreasing  tunnel  size.  For  the  15  x  15  m 
50  m  deep  case,  the  tunnel  image  is  only  19  percent  more  resistive  than  the  background 
(119  Ohm-m  maximum  resistivity  vs.  100  Ohm-m  background).  For  the  10  x  10  m  50m 
deep  case,  the  tunnel  image  is  10%  more  resistive  than  the  background  (110  Ohm-m 
maximum  resistivity  vs.  100  Ohm-m  background).  For  the  5  x  5  m  50  m  deep  case,  the 
tunnel  image  is  3%  more  resistive  than  the  background.  Hence,  the  smaller  perturbations 
from  smaller  tunnels  manifest  themselves  in  RRI  inversions  as  decreases  in  tunnel  image 
contrast  from  the  background.  These  low  contrasts  between  the  tunnel  and  the 
background  are  another  consequence  of  the  low  resolution  of  the  inversion  algorithm. 
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The  relatively  poor  resolution  exhibited  by  the  RRI  inversions  can  be  attributed  to  the 
way  in  which  RRI  computes  its  results,  and  not  to  the  inherent  resolution  limits  set  by  the 
physics  of  the  problem.  RRI  determines  a  model  structure,  m(x,z ),  for  the  conductivity  as 
a  function  of  position.  To  determine  m(x,z),  RRI  minimizes  a  functional  that  includes 
both  a  smoothness  constraint  Q  and  a  measure  of  misfit  y\  The  smoothness  constraint  Q 
favors  models  whose  conductivity  m(x,z )  has  small  second  derivatives  with  respect  to  x 
and  z.  Hence,  a  smoothly  varying  m(x,z )  with  small  second  derivatives  is  favored  over  a 
sharp  m(x,z)  with  large  second  derivatives  -  even  if  both  models  fit  the  input  data  equally 
well.  As  a  result,  RRI  tends  to  produce  smooth  images.  Although  this  approach  has  been 
useful  for  the  targets  encountered  in  geophysical  exploration,  man-made  structures  are 
smaller  and  have  sharper  contrasts,  and  the  inherent  bias  of  RRI  toward  the  smoother 
output  is  not  appropriate  for  this  type  of  target. 

It  might  be  argued  that  the  bias  toward  the  smoother  image  ensures  that  any  sharp 
features  that  are  seen  in  the  output  inversion  are  genuinely  present  in  the  input  data;  this 
point  is  correct,  but  does  not  justify  the  approach  taken  by  RRI.  The  implicit  assumption 
in  RRI’s  use  of  smoothness  is  that  smooth  conductivity  profiles  are  a  priori  more  likely 
than  rapidly  varying  profiles,  and  this  is  clearly  not  justified  by  the  characteristics  of 
man-made  structures.  Even  for  use  on  geologic  targets,  we  would  argue  that  the  tradeoff 
between  resolution  and  reliability  should  be  under  the  control  of  the  operator,  not  fixed  at 
the  most  conservative  value. 

3.2.2.  RRI  Inversion  of  Simulated  Golden  Zone  Data 

We  discuss  next  the  results  of  the  simulations  and  RRI  inversions  for  an  underground 
structure  with  characteristics  similar  to  those  of  the  Golden  Zone  survey  line  1 
experimental  data.  In  this  case  the  tunnel  depth  to  size  ratio  was  10,  similar  to  the  case 
presented  in  Figure  23.  The  tunnel  dimensions  were  2.5  x  3.0  m,  the  tunnel  depth  28  m, 
and  the  half-  space  resistivity  200  Ohm-m.  The  forward  simulation  results  for  this  case 
are  shown  in  Figure  15;  the  inversion  results  are  shown  in  Figure24.  The  inversion  mesh 
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was  similar  to  the  case  shown  in  Figure  22.  The  horizontal  sampling  at  the  surface  was 
uniform,  every  5  m,  for  the  inversion  shown  in  the  right  side  of  Figure  24,  and  was 
increased  to  every  2.5  m  between  -20  m  and  20  m  for  the  inversion  shown  show  on  the 

left  side. 


The  differences  in  the  two  inversion  images  can  be  explained  in  the  following  manner. 
When  supplied  data  every  5  meters  on  a  mesh  that  is  spaced  every  2.5  meters,  RRI 
interpolates  between  the  soundings.  However,  on  the  same  mesh,  but  with  soundings 
every  2.5  meters  in  the  region  above  the  tunnel,  interpolation  is  not  needed  above  the 
tunnel.  The  result  on  the  final  RRI  model  is  that  the  region  of  highest  conductivity 
contrast  from  the  background  is  twice  as  wide  as  the  uniformly  sampled  case.  The 
strongest  horizontal  variation  in  field  on  the  surface  occurs  in  a  neighborhood  directly 
above  the  center  of  the  tunnel.  Hence,  this  is  the  region  where  interpolation  is  the  least 
desirable.  In  short,  interpolation  of  sounding  data  near  the  center  of  the  tunnel  is  causing 
RRI  to  lose  too  much  information  about  the  slopes  of  the  field  in  this  region,  resulting  in 

a  poorer  image. 

3.2.3.  Effect  of  Surface  Layers  on  Tunnel  Detection 

Consider  next  RRI  inversions  of  tunnels  in  half-spaces  with  conducting  surface  layers. 
The  left  image  in  Figure  25  shows  the  result  of  inverting  the  forward  data  in  Figure  16. 
The  meshing  and  uniform  sampling  density,  and  the  image  intensity  scale,  are  the  same 
as  in  Figure  24.  The  tunnel,  the  intensity  peak  at  205  Ohm-m  at  a  depth  between  50  m 
and  70  m,  is  barely  visible. 
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Figure  24.  Effects  of  Sampling  Density  on  RRI  Inversion.  Inversion  of  Data 
from  Figure  15;  Sampling  Points  are  Shown  at  the  Surface. 
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Figure  25.  Effect  of  Conductive  Layer  on  RRI  Inversion  Result.  Inversions  of  Data 
from  Figure  16  With  (left)  and  Without  (right)  Tunnel. 


The  main  problem  here  is  that  the  surface  layer  is  an  order  of  magnitude  more  conductive 
than  the  background  and  therefore  introduces  a  second  contrast  scale  into  the  problem. 
On  the  other  hand,  RRI  shows  the  smaller  tunnels  as  smoothed  regions  differing  by  only 
a  few  percent  from  the  background.  Once  RRI  is  given  the  large  conductivity  contrast  due 
to  the  surface  layer,  it  does  not  perform  well  for  effects  on  a  much  smaller  conductivity 
contrast  scale.  To  highlight  the  multiple  scaling  problem,  we  also  show  the  inversion  for 
the  case  without  a  tunnel  but  with  the  same  conducting  surface  layer  and  background 
conductivity  as  in  the  previous  case.  The  result  is  shown  in  the  right  image  in  Figure  25. 
In  this  inversion,  the  image  should  be  uniform  in  x.  This  is  not  the  case.  Rather,  we  see 
banding  of  what  should  be  a  uniform  background  below  the  conducting  layer.  For  the 
case  of  uniform  half  space  with  a  tunnel  but  no  top  conducting  layer,  such  banding  does 
not  occur.  Furthermore,  we  even  see  two  perturbations  that  could  be  interpreted  as 
structures  (at  25  m,  and  -25  m,  and  both  at  a  depth  of  60  m)  even  though  there  was  no 
structure  in  this  forward  simulation. 

The  above  discussion  shows  that  the  second  conductivity  contrast  scale  introduced  by  the 
presence  of  a  conducting  surface  layer  can  mask  the  relatively  small  perturbations 
induced  by  the  tunnel.  To  correct  for  this  effect,  we  can  map  the  two-layer  problem  into 
an  equivalent  half-space.  In  the  forward  simulation,  the  equivalent  conductivity  would  be 
a  weighted  average  of  the  conductivity  of  the  top  layer,  its  thickness,  and  the  rest  of  the 
half-  space  conductivity  and  thickness.  In  the  inverse  problem,  one  could  constrain  the 
large  scale  structure  to  be  a  uniform  half-space.  Such  a  procedure  would  result  in  poor 
depth  estimates,  but  would  remove  the  large  conductivity  difference  from  the  inversion 
problem.  Once  the  tunnel  is  found,  it  is  possible  to  compute  the  background  conductivity 
more  accurately  in  order  to  estimate  the  tunnel’s  depth. 
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4.  Experiments 


A  major  part  of  this  effort  is  the  verification  of  theoretical  and  computational  studies  of 
imaging  performance  by  field  experiments.  Two  experiments  have  been  conducted  to 
date:  the  first  at  the  Superconducting  Supercollider  site  near  Waxahachee,  Texas,  and  the 
second  at  the  Golden  Zone  mine  near  Cantwell,  Alaska.  A  detailed  description  of  the 
experimental  equipment  and  procedures  is  given  in  this  section,  followed  by  an  analysis 
of  the  results  obtained  in  each  of  the  two  experiments. 

4.1.  Equipment  and  Measurement  Procedures 

A  complete  ELF/VLF  measurement  system,  including  a  transmitter,  sensors,  and  data 
recorder,  was  assembled  for  the  purpose  of  acquiring  experimental  data  under  this 
contract.  The  system  was  optimized  for  making  magnetotelluric  measurements  in  the 
VLF  frequency  range,  which  is  the  optimum  frequency  range  for  most  man-made 
underground  targets.  Each  of  the  major  components  is  described  in  detail  in  the  sections 
below. 

4.1.1.  Transmitter 

The  local  source  transmitter  was  designed  and  built  by  APTI,  and  designed  to  provide 
VLF  signals  to  supplement  natural  background  noise.  The  system  is  very  compact  when 
disassembled,  and  can  be  set  up  or  taken  down  in  approximately  30  minutes.  It  is 
composed  of  a  pair  of  triangular  vertical  loops,  supported  by  a  central  fiberglass  mast  5 
meters  in  height.  Each  of  the  two  perpendicular  loops  is  14  meters  across  at  the  base,  and 
contains  3  turns,  for  a  total  effective  area  per  loop  of  105  square  meters. 

The  two  loops  are  driven  by  a  battery  powered  two-channel  amplifier  with  a  peak  output 
voltage  of  9.0  Volts;  since  the  DC  resistance  of  each  loop  is  2.9  Ohms,  the  loop  current  at 
low  frequencies  is  3.1  Amps,  resulting  in  a  total  moment  of  325  A-m2  for  each  of  the  two 
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loops,  and  a  total  moment  of  650  A-m2.  The  inductance  of  each  loop  is  approximately  0.3 
milliHenries,  so  at  frequencies  above  1.5  kHz,  the  loop  current  decreases;  at  the  highest 
operating  frequency  of  42  kHz,  the  loop  current  has  decreased  to  0.11  Amps,  reducing  the 
total  moment  to  23  A-m  .  Since  the  power  supplied  by  the  transmitter  is  needed  most  in 
the  region  from  1  kHz  to  3  kHz  where  little  natural  noise  is  present,  the  reduction  in 
moment  with  frequency  is  not  normally  significant. 

However,  transmitter  enhancements  planned  for  the  next  year  include  the  addition  of  a 
larger  driver  amplifier  capable  of  operation  into  loads  as  small  as  1  Ohm,  and  a  meter  for 
real-time  current  monitoring.  This  will  allow  the  operator  to  adjust  the  current  level  to 
maintain  a  more  nearly  constant  moment  if  desired.  In  addition,  to  allow  operation  over 
very  resistive  ground  where  the  highest  transmitter  frequencies  are  most  important,  a 
single-turn  loop  assembly  will  be  fabricated.  This  would  reduce  the  inductance  to  0.033 
mH,  and  raise  the  frequency  at  which  moment  falloff  begins  to  14  kHz.  In  conjunction 
with  the  larger  driver  amplifier,  this  configuration  could  provide  higher  moments  at  all 
frequencies,  and  much  higher  moments  at  the  upper  end  of  the  frequency  range,  at  the 
expense  of  higher  power  consumption. 

Input  to  the  driver  amplifier  is  provided  by  a  manually  controlled  function  generator 
capable  of  producing  either  square  or  sine  wave  signals  of  adjustable  amplitude  at 
frequencies  from  10  Hz  to  80  kHz.  In  normal  operation,  sine  waves  at  the  20  frequencies 
shown  in  Table  1,  starting  at  800  Hz  and  ending  at  42  kHz,  are  used. 


56 


Table  1.  Standard  local  source  operating  frequencies. 


Frequency  Designation 
A  " 

B 
C 
D 
E 
F 
G 
H 
I 
J 
K 
L 
M 
N 
O 
P 

Q 

R 

S 

T 


Frequency  (kHz) 
0.8 
1.0 
1.2 
1.5 
2.0 
2.8 

3.2 

4.2 
5.0 

6.4 

7.5 
8.0 
10.0 
12.0 
15.0 
18.0 
20.0 
28.0 
32.0 
42.0 


Additional  transmitter  improvements  planned  for  the  next  year  include  the  ability  to  set 
transmitter  frequencies  remotely,  to  allow  the  recorder  operator  to  control  the  transmitter 
and  to  allow  computer  control  of  both  the  transmitter  and  the  recorder.  This  would 
eliminate  the  possibility  of  operator  error  in  frequency  selection  and  transmitter/recorder 
synchronization. 

4.1.2.  Sensors 

The  sensors  used  to  acquire  experimental  data  are  grouped  into  three  sensor  stations. 
Each  station  has  two  electric  field  and  three  magnetic  field  sensors,  giving  it  the 
capability  to  measure  five  of  the  six  electromagnetic  field  components  simultaneously. 
The  components  of  each  sensor  station  are  described  in  the  sections  below. 


57 


4.I.2.I.  Electric  Field  Sensors 


The  electric  field  sensors  used  to  acquire  all  experimental  data  are  models  BF-16  or  BF- 
25,  manufactured  by  ElectroMagnetic  Instruments.  The  BF-16  and  BF-25  models  are 
identical  except  for  the  length  of  the  cables  connecting  them  to  the  analog  front  end 
amplifier.  Both  models  consist  of  a  pair  of  stainless  steel  electrodes  approximately  10 
centimeters  in  length,  and  a  cable  containing  a  built-in  preamplifier.  The  measured 
frequency  response  of  the  electrode  and  preamplifier  is  shown  in  Figure  26,  and  the 
manufacturer’s  specified  noise  spectrum  is  shown  in  Figure  27. 


Figure  26.  Electric  Field  Sensor  Frequency  Response 
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Figure  27.  Specified  Electric  Field  Sensor  Noise  Spectrum 


Sensing  the  electric  field  in  a  single  direction  requires  a  pair  of  electrodes.  Both  members 
of  the  pair  are  driven  into  the  ground  a  specified  distance  apart,  and  the  voltage  difference 
is  measured;  dividing  this  voltage  by  the  known  separation  yields  the  electric  field.  The 
electric  field  sensors  for  adjacent  stations  are  positioned  so  that  the  measurements  are 
contiguous;  this  reduces  the  possibility  of  spatial  aliasing.  Normal  experimental 
procedure  is  to  space  the  stations  5  meters  apart;  although  a  wider  spacing  produces  a 
larger  signal  for  the  same  electric  field,  it  also  reduces  the  spatial  resolution  of  the 
resulting  data.  To  determine  the  effect  of  increasing  separation  on  resolution,  one  data  set 
has  been  taken  at  10  meter  spacing,  as  will  be  described  below. 

4.I.2.2.  Magnetic  Field  Sensors 

The  magnetic  field  sensors  used  are  ElectroMagnetic  Instruments  model  BF-IM,  which 
are  essentially  identical  to  the  BF-6.  The  sensors  are  induction  coil  magnetometers  5 
centimeters  in  diameter  and  80  centimeters  in  length.  They  have  a  mu-metal  core  and  an 
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integral  preamplifier  with  magnetic  feedback  for  flat  response  over  a  broad  frequency 
range.  The  specified  frequency  response  and  noise  spectrum  of  these  sensors  are  shown  in 
Figures  28  and  29.  An  additional  pair  of  BF-6-HF  sensors,  also  from  ElectroMagnetic 
Instruments,  are  used  as  a  remote  reference. 


Figure  28.  Magnetic  Field  Sensor  Frequency  Response 
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Figure  29.  Specified  Magnetic  Field  Sensor  Noise  Spectrum 


Each  sensor  station  includes  three  magnetic  field  sensors,  allowing  simultaneous 
measurements  of  the  magnetic  field  in  all  three  axes.  To  date,  all  experimental  data  has 
included  only  the  horizontal  electric  and  magnetic  field  components,  principally  because 
of  the  additional  setup  time  required  for  vertical  magnetic  measurements.  Except  in  still 
conditions,  measurement  of  the  vertical  magnetic  field  requires  burying  the  vertical 
magnetic  sensor  to  prevent  wind-induced  vibration.  This  slows  the  process  of  station 
setup  significantly.  In  addition,  simulation  results  suggest  that  the  addition  of  a  vertical 
magnetic  field  measurement  is  most  helpful  on  highly  conductive  targets;  since  all  of  the 
targets  examined  to  date  have  been  resistive,  it  is  expected  that  the  absence  of  a  vertical 
magnetic  measurement  has  not  significantly  reduced  imaging  performance.  Experiments 
on  conductive  targets,  including  vertical  magnetic  field  measurements,  are  planned  for  the 
next  year. 
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4.I.2.3.  Analog  Front  End  Amplifier 


Both  the  electric  and  magnetic  field  sensor  outputs  are  passed  through  the  analog  front 
end  amplifier  before  being  sent  to  the  data  recorder.  The  principal  purpose  of  the  analog 
front  end,  also  produced  by  Electromagnetic  Instruments,  is  to  provide  additional 
amplification  and  filtering  of  power  line  signals  and  their  harmonics.  The  electric  field 
channels  in  the  analog  front  end  have  an  adjustable  gain  from  10  to  3200,  the  magnetic 
field  channels  from  1  to  800.  All  five  channels  also  contain  switchable  highpass,  lowpass, 
and  power  line  filters;  in  normal  operation,  the  filters  are  set  to  provide  flat  response  and 
filtering  of  60,  120  and  180  Hz  power  line  signals.  The  measured  response  of  an  electric 
field  channel  in  this  configuration  is  shown  in  Figure  30. 


Figure  30.  Analog  Front  End  Frequency  Response 
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4.I.2.4.  Data  Recorder 


Signals  from  each  of  the  three  stations,  as  well  as  the  two  remote  magnetic  reference 
channels,  are  recorded  by  a  Storeplex  Delta  digital  data  recorder  manufactured  by  Racal 
Recorders.  By  recording  all  of  the  time  series  data,  rather  than  just  spectra,  it  is  possible 
to  examine  the  recorded  data  in  detail  and  process  it  with  several  different  techniques  to 
compare  the  results.  The  recorder  is  equipped  with  20  analog  input  channels,  each  with  a 
16-bit  sigma-delta  analog-to-digital  converter  capable  of  operating  at  a  sampling  rate  of 
up  to  100  kHz.  The  sigma-delta  converters  used  in  this  system  contain  internal  digital 
decimators  that  reduce  an  internal  6.4  megahertz  low  resolution  output  sample  stream  to  a 
higher  resolution  100  kHz  output  sample  rate,  while  applying  a  linear  phase  digital 
antialiasing  filter  with  ±  0.1  dB  ripple  from  0  to  45.5  kHz  and  greater  than  90  dB  of 
stopband  attenuation.  A  significant  advantage  of  sigma-delta  converters  is  that  the  digital 
antialiasing  filter  eliminates  the  need  for  all  but  a  simple  second  order  analog  filter  at  the 
recorder  input. 

The  recorder  is  capable  of  recording  14  channels  of  100  kHz  sampled  data  (two  electric 
and  two  magnetic  channels  for  each  of  the  three  sensor  stations,  and  two  remote  reference 
channels)  for  up  to  60  minutes  on  a  single  ST-160  S-VHS  video  tape.  When  three 
additional  vertical  channels  are  added,  the  maximum  recording  time  is  reduced  to  30 
minutes  at  a  100  kHz  sampling  rate.  The  recorder  can  also  sample  at  6.25,  12.5,  25,  and 
50  kHz  when  high  frequency  data  is  not  required,  with  corresponding  increases  in  the 
recording  time.  An  additional  advantage  of  the  sigma-delta  converters  is  that  the  digital 
antialiasing  filters  automatically  adapt  to  changes  in  the  sampling  rate. 

The  dynamic  range  of  a  16-bit  sampling  system  is  somewhat  smaller  than  the  specified 
dynamic  range  of  the  sensors,  so  care  must  be  taken  to  match  the  analog  front  end  gains 
to  the  signal  amplitude  so  that  quantization  noise  does  not  contribute  to  the  noise  floor  of 
the  overall  data  acquisition  system.  The  ideal  noise  floor  of  a  perfect  16-bit  data 
acquisition  system  and  the  measured  noise  floor  of  the  recorder  are  shown  in  Figure  31. 


63 


The  ideal  system  has  a  dynamic  range  of  96  dB;  the  recorder  delivers  approximately  90 
dB  of  dynamic  range,  a  very  good  performance  for  a  system  sampling  at  the  relatively 
high  rate  of  100  kHz. 


Figure  31.  Ideal  and  Measured  Recorder  Noise  Spectra. 


A  comparison  of  the  noise  contribution  of  the  recorder  and  the  intrinsic  noise  of  a 
magnetic  field  sensor  with  a  typical  field  spectrum  is  shown  in  Figure  32.  The  signal 
spectrum  shown  is  recorded  from  one  of  the  remote  reference  sensors,  which  do  not  pass 
through  the  analog  front  end  and  have  no  filtering  of  power  line  harmonics.  This 
represents  the  worst  case  with  respect  to  the  dynamic  range  demands  placed  upon  the  data 
recorder.  Note  that  in  this  case,  both  the  recorder’s  quantization  noise  and  the  intrinsic 
noise  of  the  sensor  are  well  below  the  level  of  the  recorded  noise  spectrum,  indicating 
that  the  recorded  signal  is  an  accurate  representation  of  the  actual  magnetic  field  level. 
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Figure  32.  Comparison  of  Recorder  and  Sensor  Noise  with  Typical  Measured  Signal 

In  examining  the  recorder  and  sensor  noise  spectra  in  Figure  32,  we  can  see  that  the 
quantization  noise  of  the  recorder  is  somewhat  higher  than  the  intrinsic  noise  of  the 
sensor  in  the  region  from  1  to  10  kHz;  this  would  be  the  limiting  factor  in  this  frequency 
region  if  the  natural  magnetic  signal  level  were  very  low.  The  sensor’s  internal  noise  is 
fixed  by  the  design  of  the  sensor,  and  cannot  be  adjusted,  but  if  the  signal  level  is  low,  it 
is  possible  to  adjust  the  recorder  gain  so  that  the  quantization  noise  is  reduced.  This 
procedure  has  been  followed  whenever  possible. 

The  data  recorder  also  automatically  records  date  and  time  information  accurate  to  0.01 
seconds  for  all  recorded  data.  Before  each  recording  session,  the  recorder  clock  is 
manually  synchronized  to  the  UTC  time  obtained  from  a  handheld  GPS  receiver,  ensuring 
that  the  recorded  time  is  accurate  to  approximately  ±0.1  seconds.  If  higher  accuracy  is 
needed  in  the  future,  the  recorder  can  also  accept  an  IRIG-B  timecode  input  that  can  be 
produced  by  a  portable  GPS  time  standard,  allowing  the  recorded  time  to  be  matched  to 
UTC  to  within  the  10  millisecond  resolution  of  the  recorder  clock. 
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4.I.2.5.  Data  Transfer 


The  data  recorder  includes  a  SCSI-2  interface  to  allow  direct  transfer  of  the  recorded  data 
to  a  computer  for  analysis.  The  data  are  transferred  in  16-bit  digital  form  under  the 
control  of  data  transfer  software  provided  by  the  manufacturer.  During  the  course  of 
processing  experimental  data,  it  was  discovered  that  the  data  received  by  the  computer 
were  corrupted  by  occasional  errors.  These  errors,  although  infrequent,  were  sufficient  to 
corrupt  the  spectra  and  make  them  unusable  for  imaging.  In  order  to  obtain  uncorrupted 
spectra,  it  was  necessary  to  repeat  the  data  transfer  process. 

In  examining  this  problem,  it  was  determined  that  the  fault  lay  with  the  data  transfer 
software.  The  errors  appearing  in  the  outputs  were  not  flaws  in  the  data  recording  tape, 
because  they  were  not  replicated  when  a  second  transfer  of  the  same  section  of  tape  was 
made.  To  rule  out  a  hardware  flaw,  an  identical  data  recorder  was  obtained  and  tested;  the 
results  were  similar  using  both  recorders.  Different  computers  and  SCSI  controllers  were 
also  tested,  with  similar  results.  The  only  common  element  was  the  manufacturer- 
provided  data  transfer  software.  A  preliminary  version  of  the  next  software  revision  was 
obtained  from  the  manufacturer;  this  version  reduced  but  did  not  eliminate  the  data 
transfer  errors. 

Although  the  errors  could  be  corrected  by  repeating  the  data  transfer  process,  this 
problem  made  it  necessary  to  manually  examine  the  output  of  each  data  transfer  to 
determine  which  were  correct  and  which  should  be  repeated.  Because  thousands  of 
spectra  are  computed  for  each  line  of  sensor  data,  the  process  of  manual  examination 
greatly  slowed  the  analysis  of  experimental  results. 

Since  this  problem  has  had  a  serious  impact  on  the  efficiency  of  the  analysis  process,  a 
permanent  solution  has  been  developed.  A  high-speed  digital  interface  that  bypasses  the 
SCSI  subsystem  and  does  not  require  use  of  the  manufacturer’s  software  has  been  ordered 
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and  is  currently  under  test.  This  hardware  solution  should  solve  the  transfer  error  problem 
and  remove  the  need  for  manual  examination  of  spectra.  The  new  hardware  is  expected  to 
be  delivered  before  the  start  of  the  next  field  experiment  and  will  be  utilized  in  all  future 
processing. 

4.1.3.  Standard  Procedures 

The  normal  data  acquisition  procedure  is  to  place  three  stations  at  adjacent  locations  5 
meters  apart,  and  record  electric  and  magnetic  field  levels  at  these  three  stations 
simultaneously.  A  remote  reference  station,  composed  of  two  perpendicular  magnetic 
field  sensors,  is  also  recorded  at  the  same  time  as  the  station  signals.  The  remote 
reference  is  typically  placed  70-100  meters  from  the  sensors,  and  is  connected  to  the  data 
recorder  by  a  shielded  coaxial  cable. 

To  reduce  the  amount  of  wind-induced  noise  due  to  motion  of  the  magnetic  sensors,  these 
sensors  are  stabilized  with  sandbags,  or,  under  windy  conditions,  buried  in  trenches. 
Coupling  of  wind-induced  vibration  to  the  sensor  outputs  is  less  severe  at  the  VLF 
frequencies  of  interest  here  than  at  ELF  frequencies;  nevertheless,  it  can  be  the  principal 
contributor  to  the  noise  background  if  care  is  not  taken  to  stabilize  the  sensors. 

If  the  ground  conditions  are  dry,  the  electrodes  used  for  electric  field  measurements  are 
wet  with  water  to  reduce  the  contact  resistance.  As  in  the  case  of  wind-induced  noise  for 
the  magnetic  field  sensors,  noise  caused  by  the  contact  resistance  of  the  electrodes  is  less 
of  a  problem  at  VLF  than  at  lower  frequencies.  In  both  of  the  experiments  described 
below,  the  ground  surface  was  sufficiently  wet  that  additional  moisture  was  not  needed. 

Although  the  sensor  station  configurations  are  identical,  the  data  recording  procedure 
varies  depending  on  whether  a  local  source  or  an  ionospheric  source  is  used.  The  standard 
procedures  for  each  of  these  are  described  below. 
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4.1.3.L  Local  Source 


When  the  transmitter  described  above  is  used  as  a  local  source,  twenty  minutes  of  data 
are  recorded  at  each  station.  One  minute  of  data  is  recorded  at  each  of  the  transmitter 
frequencies  shown  in  Table  1.  To  allow  accurate  identification  of  the  frequencies  using 
the  time  data  recorded  automatically  by  the  recorder,  the  transmitter  frequency  schedule  is 
synchronized  to  the  recorder  clock,  so  that  all  frequency  transitions  occur  between  :50 
and  :00  seconds  on  the  recorder  clock.  This  permits  the  transmitter  frequencies 
corresponding  to  each  minute  of  recorded  data  to  be  easily  identified,  and  ensures  that  at 
least  50  seconds  of  usable  data  are  obtained  at  each  frequency. 

The  range  of  800  Hz  to  42  kHz  covered  by  the  standard  local  source  frequency  schedule 
is  appropriate  for  a  wide  range  of  ground  conductivities,  and  it  has  not  been  found 
necessary  to  modify  this  schedule  for  any  of  the  experiments  performed  to  date.  The 
remainder  of  the  measurement  system  limits  the  useful  frequency  range  to  below  50  kHz; 
this  limit  is  set  both  by  the  rising  noise  level  of  the  magnetic  field  sensors  and  the  45.5 
kHz  cutoff  frequency  of  the  recorder.  As  described  in  the  transmitter  section  above,  some 
source  modifications  are  possible  to  increase  the  transmitted  power  at  high  frequencies; 
this  would  increase  the  available  power  at  frequencies  up  to  the  45.5  kHz  recorder  cutoff. 

At  frequencies  below  1.5  kHz,  the  output  of  the  transmitter  is  constant  with  frequency 
down  to  the  lower  frequency  limit  of  10  Hz.  The  major  difficulty  in  employing  the  local 
source  at  low  frequencies  is  placing  it  sufficiently  far  from  the  measurement  site  so  that 
its  signal  approximates  that  of  a  distant  source.  In  order  to  meet  this  requirement,  the 
transmitter  must  be  positioned  at  least  3  skin  depths  (in  the  ground  at  the  experiment 
location)  away  from  the  measurement  location.  As  the  conductivity  of  the  ground  drops, 
this  distance  becomes  larger,  and  since  the  field  from  the  transmitter  falls  off  as  1/r3,  it 
becomes  difficult  to  obtain  sufficient  signal  from  the  transmitter  at  the  distances  that  must 
be  employed  in  resistive  grounds. 


68 


Although  these  factors  could  limit  the  utility  of  the  local  source  in  regions  of  very  high  or 
very  low  ground  conductivity,  little  difficulty  has  been  experienced  with  these  effects  in 
the  experiments  conducted  to  date.  In  some  cases,  the  experimental  data  show  that  the  3 
skin  depth  minimum  distance  has  not  been  met  for  the  lower  end  of  the  transmitter 
frequency  range.  This  is  not  so  much  a  problem  of  transmitter  design  as  of  logistics;  the 
transmitter  location  is  also  affected  by  practical  considerations  such  as  ease  of  operator 
access.  The  modifications  to  increase  transmitter  power  described  in  the  transmitter 
section  above  should  give  increased  flexibility  in  selecting  the  transmitter  location. 

4.I.3.2.  Ionospheric  Source 

When  an  ionospheric  source  is  used,  the  signals  from  the  source  are  typically  much 
weaker  than  those  from  a  local  transmitter.  To  accommodate  the  reduced  signal  level,  the 
amount  of  time  spent  at  each  frequency  is  doubled  to  two  minutes,  and  a  total  of  forty 
minutes  of  data  are  recorded  at  each  sensor  location.  Synchronization  of  the  recorder 
clock  with  the  transmitted  signal  is  made  possible  by  setting  the  recorder  clock  using  a 
GPS  derived  time  signal  as  described  above.  The  ionospheric  heater’s  transmission 
schedule  is  also  accurately  synchronized  to  GPS  time,  so  the  location  of  the  frequency 
transitions  of  the  ionospheric  source  can  be  identified  with  an  error  of  much  less  than  a 
second.  As  with  the  local  source,  this  allows  easy  identification  of  the  transmitter 
frequency  using  time  data  recorded  along  with  the  measured  signals. 

The  frequency  range  used  for  the  ionospheric  source  is  selected  for  each  experimental  site 
on  the  basis  of  conductivity  estimates  derived  from  local  source  measurements.  The 
ionospheric  source  does  not  suffer  from  the  low  frequency  limitations  that  the  local 
transmitter  does,  since  it  is  always  at  a  large  distance  from  the  measurement  site.  At 
higher  frequencies,  an  ionospheric  source  has  approximately  constant  signal  output  (for 
fixed  ionospheric  conditions)  up  to  the  point  at  which  it  can  no  longer  raise  the 
ionospheric  electron  temperature  fast  enough  to  keep  pace  with  the  modulating 
waveform.  This  limit  is  determined  by  the  RF  power  density  achieved  at  the  heating 
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altitude  by  the  ionospheric  source,  and  depends  on  the  heating  frequency,  gain  and 
transmitter  power  of  the  heating  facility.  For  the  HAARP  DP  facility  at  330  kW  RF 
power,  it  is  approximately  5  kHz,  for  the  HIP  AS  facility,  approximately  8  kHz,  and  for 
the  completed  HAARP  facility  at  3.6  MW,  it  is  estimated  to  be  20  -  40  kHz. 

4.I.3.3.  Surveying 

To  provide  reliable  data  on  the  location  of  the  sensor  stations  relative  to  the  target,  the 
locations  are  surveyed  using  a  total  station  instrument.  Survey  data  are  adjusted  to 
produce  accurate  position  estimates  and  ensure  consistency  of  the  measured  angles  and 
distances;  the  typical  residual  error  (the  estimated  standard  deviation  of  the  adjusted 
points)  is  less  than  0.05  meters.  In  many  cases,  the  target  itself  is  not  accessible,  so  data 
must  be  obtained  on  the  location  of  the  target  relative  to  surface  landmarks  that  can  be 
surveyed.  In  this  case,  the  accuracy  of  location  estimates  also  depends  on  the  accuracy  of 
the  target  location  data. 

4.2.  SSC  Experiment 

Two  trips  to  perform  experiments  at  the  Superconducting  Supercollider  site  near 
Waxahatchee,  Texas  were  made  on  March  10-14  and  April  7-11,  1997.  The  original 
intent  was  to  acquire  all  the  necessary  data  in  one  trip;  however,  because  of  poor  weather, 
very  little  data  was  acquired  during  the  first  visit  in  March.  On  the  second  trip,  it  was 
decided  to  repeat  the  acquisition  of  the  data  acquired  on  the  first  trip  to  ensure  that  all  the 
data  were  gathered  under  similar  conditions. 
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4.2.1.  Description  of  Experiment 


The  target  for  this  experiment  was  the  main  ring  of  the  Superconducting  Supercollider  at 
a  location  just  north  of  the  building  complex  at  station  N-15.  At  this  location,  the  main 
ring  is  5.6  meters  in  diameter  and  approximately  65  meters  beneath  the  surface.  This 
section  of  the  tunnel  is  lined  with  reinforced  concrete  in  prefabricated  sections.  Because 
pumping  of  water  from  the  tunnel  ceased  more  than  two  years  before  this  experiment,  it  is 
assumed  that  the  tunnel  is  partially  or  completely  filled  with  groundwater. 

Data  were  taken  along  two  lines  of  length  105  meters,  as  shown  in  Figure  33.  Line  1 
crossed  the  tunnel  approximately  perpendicular  to  the  tunnel  centerline;  it  was  surveyed 
both  at  the  usual  5  meter  station  spacing  and  at  a  spacing  of  10  meters  to  determine  the 
effects  of  station  spacing  on  resolution.  The  second  line  crossed  the  tunnel  at 
approximately  a  45°  angle  to  the  tunnel  centerline  and  was  taken  to  examine  the  effects  of 
crossing  angle  on  resolution  of  two-dimensional  targets. 

The  line  locations  were  surveyed  with  respect  to  the  intersections  of  Loop  Road  with 
Hoyt  Road,  and  Hoyt  Road  with  FM  1446,  as  well  as  the  comers  of  the  MTL  building 
and  the  surface  mounds  associated  with  the  personnel,  utility,  and  magnet  delivery  shafts. 
The  end  points  of  a  USGS  survey  conducted  during  the  first  visit  were  also  surveyed  and 
are  shown  on  Figure  33.  These  locations  were  correlated  with  data  obtained  later  to 
accurately  locate  the  tunnel;  it  is  estimated  that  the  horizontal  location  of  the  tunnel  is 
known  to  within  1  meter,  and  the  depth  to  within  5  meters. 
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SSCL  N-15  SITE  MAP 
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Figure  33.  Site  Map  of  the  SSC  Experiment. 
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4.2.2.  Analysis  of  Results 


As  the  data  from  the  survey  were  analyzed,  it  became  apparent  that  one  of  the  three 
sensor  stations  exhibited  a  systematic  gain  difference  that  was  consistent  throughout  the 
experiment.  Although  examination  of  the  time  series  data  showed  that  the  data  were 
valid,  the  response  of  the  station  was  different  from  the  other  two  stations  by  a  significant 
amount.  The  size  of  the  difference  and  the  fact  that  it  appeared  at  every  location  where 
that  sensor  station  was  employed  ruled  out  the  possibility  that  the  variation  resulted  from 
genuine  conductivity  variations  in  the  ground. 

In  order  to  correct  for  this  error,  the  data  from  the  station  in  question  were  normalized  so 
that  their  average  apparent  resistivity  along  the  entire  survey  line  was  similar  to  the  other 
two  stations.  This  reduced  the  variation  sufficiently  to  allow  the  data  to  be  analyzed,  but 
did  not  completely  eliminate  the  problem.  Extensive  calibration  tests  of  the  sensors  did 
not  detect  any  malfunction  in  the  sensors  or  the  analog  front  end;  it  appears  that  the 
difference  must  have  resulted  from  a  systematic  error  in  setting  up  the  station  or  in 
programming  its  gain  at  each  of  the  sensor  positions. 

Shown  in  Figure  34  is  an  inversion  result  obtained  using  the  5-meter  spacing  data  from 
survey  line  1.  The  true  location  of  the  tunnel  is  between  stations  2  and  3,  at  a  depth  of 
approximately  65  meters.  This  result  was  a  disappointment  for  a  number  of  reasons,  the 
principal  one  of  course  being  that  no  identifiable  response  at  the  known  location  of  the 
tunnel  was  observed.  The  residual  effects  of  the  partially  corrected  gain  variation  can  be 
seen  as  vertical  striations.  The  large  resistive  anomaly  seen  at  stations  5-11  is  similar  to 
the  expected  signature  of  the  tunnel,  but  this  would  require  an  error  in  position  location  of 
approximately  25  meters,  which  is  extremely  unlikely.  It  is  more  likely  that  this  resistive 
anomaly  is  not  the  tunnel  itself,  but  perhaps  caused  by  a  variation  in  groundwater  flow 
due  to  the  presence  of  the  tunnel. 
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Figure  34.  Resistivity  Image  for  SSC  Experiment  using 
Parameters  k  =  .77,  Xv  =  .01,  Xh  =  50;  %(77)  =  24.0. 
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The  parameters  which  affect  the  performance  of  the  imaging  algorithm  have  been 
described  in  Section  2.2.  The  algorithm  parameters  used  for  this  result  were  a  minimum 
coherence  threshold  of  K=.ll,  a  spectral  regularization  parameter  of  A =.01,  and 
horizontal  Lagrange  multiplier  of  A  =50.  This  coherence  threshold  delivered  a  percent 
over  threshold  of  %(.77)=24.  This  is  an  unusually  small  k  compared  to  what  is  usually 
employed,  but  was  chosen  to  mitigate  the  effects  of  the  gain  variation  discussed  above. 
Both  the  spectral  regularization  parameter  A  and  the  horizontal  parameter  \  are 
significantly  smaller  than  usual. 

During  the  analysis  of  the  Golden  Zone  data  described  below,  several  modifications  to  the 
processing  and  inversion  procedures  were  implemented.  These  changes  resulted  in 
significantly  better  performance  on  those  data,  and  data  from  the  SSC  survey  will  be 
reprocessed  using  the  modified  procedures  to  determine  whether  better  results  can  be 
achieved.  Processing  of  the  remaining  lines  from  the  SSC  survey  will  also  be  completed 
using  the  updated  procedures. 

4.3.  Golden  Zone  Experiment 

A  full  experimental  campaign,  including  use  of  the  HAARP  ionospheric  source,  was 
conducted  at  the  Golden  Zone  mine  near  Cantwell,  Alaska,  on  August  4-18,  1997.  The 
details  of  this  experiment  and  an  analysis  of  the  results  are  given  in  the  sections  below. 

4.3.1.  Description  of  Experiment 

The  Golden  Zone  is  a  gold  mine  whose  initial  workings  date  from  the  early  part  of  this 
century;  it  was  last  extensively  mined  in  the  1930s.  Four  levels  of  workings  exist,  at 
depths  of  approximately  100,  200,  325,  and  500  feet  from  the  top  of  the  ore  body.  The 
depth  below  ground  level  varies  depending  on  the  location. 

Of  the  four  levels,  the  325  foot  level  was  the  most  recent,  having  been  opened  in  1983, 
and  was  the  only  level  accessible  for  inspection.  The  tunnel  is  approximately  3.0  meters 
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in  diameter,  and  contains  two  compressed  air  lines  approximately  10  cm  in  diameter;  over 
most  of  its  length  the  roof  of  the  tunnel  has  been  lined  with  wire  mesh.  A  typical  interior 
view  of  the  325  foot  level  is  shown  in  Figure  35. 

A  plan  view  of  the  Golden  Zone  experiment  setup  is  shown  in  Figure  36.  Data  were  taken 
along  three  lines;  all  three  lines  passed  over  the  325  foot  level  tunnel,  at  depths  of  29 
meters  for  survey  line  1,  36  meters  for  survey  line  2,  and  48  meters  for  survey  line  3.  All 
three  lines  were  surveyed  using  the  local  transmitter  at  the  position  shown,  following  the 
standard  procedures  described  above.  In  addition,  survey  line  1  was  surveyed  on  seven 
consecutive  nights,  12  August  -  18  August  universal  time  (11  August  -  17  August  local 
time),  using  the  HAARP  facility  as  a  transmitter. 

The  HAARP  facility  was  248  kilometers  from  the  experimental  site  at  a  bearing  of  1 10° 
true.  On  each  night,  the  facility  transmitted  in  X  polarization  at  an  RF  power  level  of 
approximately  300-320  kW.  The  HF  carrier  was  amplitude  modulated  with  a  modulation 
level  of  90%  by  a  sinusoidal  signal  at  the  desired  VLF  frequency.  The  modulation 
frequency  schedule  is  shown  in  Table  2  below. 
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Figure  35.  Golden  Zone  Mine  Interior 
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Figure  36.  Golden  Zone  Mine  Plan  View. 


Table  2.  HAARP  Modulation  Frequency  Schedule 
Minutes  After  Start _ Modulation  Frequency  (kHz) 


00  -  :02 

1.0 

02  -  :04 

1.1 

04  -  :06 

1.3 

06  -  :08 

1.4 

08  -  :10 

1.6 

10  -  :12 

1.8 

12  -  :14 

2.0 

14  -  :16 

2.3 

16  - :  18 

2.6 

1 8  —  :20 

3.0 

20  -  :22 

3.3 

22  -  :24 

3.8 

24  -  :26 

4.3 

26  -  :28 

4.8 

28  -  :30 

5.5 

30  -  :32 

6.1 

32  -  :34 

7.0 

34  -  :36 

7.8 

36  -  ;38 

8.9 

38  -  :40 

10.0 

40  -  :00 

OFF 

This  schedule  was  repeated  five  times  each  night  at  one  hour  intervals.  During  each  hour, 
data  was  collected  at  three  locations  on  survey  line  1,  allowing  all  15  stations  on  that 
survey  line  to  be  covered  each  night.  The  starting  time  for  HAARP  operation  varied  over 
the  seven  nights  on  which  data  were  collected;  the  times  for  each  night  are  shown  in 
Table  3  below. 

Table  3.  HAARP  Starting  Times 


Night  1 

20:00  11  Aug  1997  ADT 

04:00  12  Aug  1997  UT 

Night  2 

20:30  12  Aug  1997  ADT 

04:30  13  Aug  1997  UT 

Night  3 

20:30  13  Aug  1997  ADT 

04:30  14  Aug  1997  UT 

Night  4 

22:30  14  Aug  1997  ADT 

06:30  15  Aug  1997  UT 

Night  5 

22:00  15  Aug  1997  ADT 

06:00  16  Aug  1997  UT 

Night  6 

22:00  16  Aug  1997  ADT 

06:00  17  Aug  1997  UT 

Night  7 

22:00  17  Aug  1997  ADT 

06:00  18  Aug  1997  UT 
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Table  4.  HAARP  Carrier  Frequencies  and  RF  Powers 
Time  (UT) _ Carrier  Frequency  (MHz)  RF  Power  (kW) 


12  Aug  1997  (Night  1) 


04:00-04:40 

3.200 

340 

05:00-05:40 

3.200 

340 

06:00  -  06:40 

3.200 

330 

07:00  -  07:40 

3.200 

320 

08:00-08:40 

3.200 

320 

13  Aug  1997  (Night  2) 
04:30-05:10 

4.438 

300 

05:30-06:10 

4.438 

320 

06:30-07:10 

4.438 

320 

07:30-08:10 

4.438 

300 

08:30-09:10 

4.438 

300 

14  Aug  1997  (Night  3) 
04:30-05:10 

4.438 

300 

05:30-06:10 

4.438 

300 

06:30-07:10 

4.438 

300 

07:30-08:10 

4.438 

300 

08:30-09:10 

4.438 

300 

15  Aug  1997  (Night  4) 
06:30-07:10 

4.438 

280 

07:30-08:10 

4.438 

300 

08:30-09:10 

4.438 

300 

09:30-10:10 

4.438 

300 

10:30-11:10 

4.438 

300 

16  Aug  1997  (Night  5) 
06:00-06:40 

4.438 

320 

07:00  -  07:40 

4.438 

320 

08:00-08:40 

4.438 

320 

09:00-09:40 

4.438 

320 

10:00-10:40 

4.438 

300 

17  Aug  1997  (Night  6) 
06:00-06:40 

4.438 

320 

07:00  -  07:40 

4.438 

320 

08:00-08:40 

4.438 

320 

09:00-09:40 

4.438 

320 

10:00-10:40 

4.438 

320 

18  Aug  1997  (Night  7) 
06:00-06:40 

4.438 

320 

07:00-07:40 

4.438 

320 

08:00-08:40 

4.438 

320 

09:00-09:40 

4.438 

320 

10:00-10:40 

4.438 

320 
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There  was  a  slight  variation  in  the  transmitted  power  level  as  the  experiment  progressed 
due  to  changes  in  the  number  of  transmitters  operational.  In  addition,  the  carrier 
frequency  was  different  on  the  first  night  of  operation.  The  complete  schedule  of  carrier 
frequencies  and  RF  power  levels  for  each  hour  of  transmission  is  shown  in  Table  4. 

On  16  August  (Night  5),  difficulties  with  the  transmitter  control  system  at  the  HAARP 
facility  caused  some  of  the  desired  frequencies  not  to  be  transmitted.  The  scheduled  start 
times  for  the  five  hours  of  operation  were  06:00,  07:00,  08:00,  09:00,  and  10:00  UT;  the 
actual  start  times  were  06:05, 07:25;  08:00,  09:00,  and  10:10  UT.  For  the  three  hours  that 
the  actual  start  times  were  missed,  the  frequencies  that  would  have  been  transmitted  were 
skipped,  so  that  the  remaining  frequencies  were  transmitted  according  to  the  schedule  in 
Table  2. 

During  this  experiment,  the  HAARP  facility  was  operating  at  a  very  low  power  level; 
consequently,  the  signal  received  at  the  experiment  site  was  relatively  small.  Because  of 
difficulties  encountered  in  transferring  data  from  the  data  recorder,  as  discussed  in 
Section  4. 1.2.5,  not  all  of  the  HAARP  data  recorded  have  been  analyzed  at  this  time. 
However,  sections  of  the  data  spaced  approximately  one  hour  apart  have  been  analyzed, 
and  it  is  possible  to  give  general  estimates  of  the  signal  levels  obtained. 

The  minimum  signal  level  useful  for  imaging  at  the  Golden  Zone  site  would  have  been 
approximately  0.05  pT  RMS;  although  the  HAARP  signal  was  often  easily  detectable  in 
the  recorded  data,  there  were  few  instances  when  the  signal  reached  0.05  pT.  Diagnostics 
available  at  the  HAARP  site  indicated  that  ionospheric  conditions  for  VLF  generation 
ranged  from  moderate  to  poor  over  the  period  of  the  experiment.  Unfortunately,  no  real 
time  data  from  local  receivers  were  available  to  the  HAARP  operator  to  allow  selection 
of  the  optimum  RF  frequency,  so  it  is  impossible  to  know  whether  higher  signal  levels 
could  have  been  achieved  by  heating  at  different  frequencies. 
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The  maximum  signal  obtained  from  the  HAARP  transmitter  in  the  data  records  analyzed 
to  date  was  approximately  0.1  pT  RMS.  This  signal  was  obtained  for  a  brief  period  at 
approximately  06:48  UT  on  13  August,  while  HAARP  was  transmitting  at  a  modulation 
frequency  of  3  kHz.  A  field  spectral  density  plot  of  this  signal,  computed  from  data  taken 
at  06:48:00  to  06:50:00  UT,  is  shown  in  Figure  37.  The  origin  of  this  signal  is  confirmed 
by  the  fact  that  the  previous  and  following  two  minute  data  records  show  similar  strong 
signals  at  2.6  kHz  and  3.3  kHz,  respectively,  the  corresponding  HAARP  modulation 
frequencies  for  those  times. 


Frequency  (Hz) 


Figure  37.  Field  Spectrum  Showing  HAARP  Signal  at  3.0  kHz 

In  general,  it  appears  that  the  signals  obtained  from  HAARP  were  not  sufficiently  strong 
to  be  useful  for  imaging  over  most  of  the  experimental  period.  For  the  facility  at  a  power 
level  of  320  kW  to  produce  signals  strong  enough  to  be  useful  at  the  250  kilometer 
distance  of  the  experimental  site,  good  to  excellent  ionospheric  conditions  would  be 
required.  The  brief  periods  of  strong  signal  obtained  indicate  that  these  conditions  do 
occur,  but  it  appears  that  during  the  late  summer  period  they  do  not  occur  frequently.  It  is 
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unfortunate  that  the  facility  operators  did  not  have  real  time  diagnostics  of  the  VLF  signal 
level  available  to  allow  them  to  select  an  optimum  heating  frequency. 

The  signal  level  obtained  from  the  HAARP  facility  operating  at  a  heating  frequency  of 
4.438  MHz  appears  to  have  smaller  amplitude  variation  from  night-to-night  than  has  been 
previously  observed  in  experiments  using  the  HIPAS  facility  at  a  heating  frequency  of 
2.85  MHz.  Although  this  observation  is  based  on  small  sets  of  data  taken  at  different 
times  and  under  different  conditions,  it  is  in  accord  with  some  previous  theoretical 
analyses  of  heating.  An  experiment  to  test  the  relative  efficiency  and  stability  of  heating 
at  different  RF  frequencies  could  easily  be  conducted  using  the  HAARP  facility  and  a 
single  set  of  VLF  receiving  coils.  Because  of  the  importance  of  signal  level  to  the 
imaging  mission,  plans  are  underway  to  conduct  such  an  experiment  during  the  next  year; 
this  would  provide  quantitative  data  to  determine  how  to  optimize  the  ionospheric 
generation  of  VLF  signals. 

The  HAARP  facility  is  in  the  process  of  being  upgraded  to  a  RF  power  level  of  960  kW; 
it  is  currently  expected  that  the  facility  will  be  available  for  experiments  at  this  power 
level  in  the  Fall  of  1998.  The  threefold  increase  in  RF  power  level  is  expected  to  result  in 
approximately  a  threefold  increase  in  field  levels  (a  factor  of  nine  increase  in  radiated 
VLF  power).  Scaling  the  data  observed  at  the  Golden  Zone  to  account  for  the  increased 
power  level,  it  appears  that  useful  signals  would  have  been  observed  as  much  as  50%  of 
the  time  had  the  960  kW  facility  been  available  for  this  experiment. 

4.3.2.  Analysis  of  Results 

This  section  contains  an  analysis  of  data  from  Golden  Zone  survey  line  1,  taken  using  the 
local  transmitter.  Data  from  local  transmitter  measurements  on  survey  lines  2  and  3,  as 
well  as  from  the  seven  nights  of  data  collection  using  the  HAARP  source,  are  currently 
being  analyzed. 
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One  effect  noted  in  the  initial  analysis  of  the  survey  line  1  data  was  that  at  the  lowest 
source  frequencies,  the  transmitter  was  too  close  to  provide  an  approximation  to  a  distant 
plane  wave.  Figure  36  shows  that  line  1  was  the  closest  to  the  transmitter  of  all  the  lines; 
source  proximity  effects  should  be  much  smaller  at  the  other  locations.  Analysis  of  the 
data  showed  that  at  line  1,  the  transmitter  signal  was  usable  only  at  frequencies  of  6.4 
kHz  and  above.  All  data  at  frequencies  below  this  were  obtained  using  natural 
background  noise  recorded  concurrently  with  the  higher  transmitter  frequencies  to  ensure 
that  the  plane  wave  assumption  inherent  in  the  inversion  was  obeyed. 

In  order  to  illustrate  the  factors  affecting  the  output  image,  we  will  first  examine  several 
different  images  derived  from  the  same  input  data  set.  We  will  begin  by  examining 
images  from  the  survey  line  1  local  source  data  as  key  parameters  of  the  algorithm  are 
adjusted,  to  provide  insight  into  the  effects  of  these  parameters.  All  detection,  estimation 
and  inversion  algorithms  have  similar  thresholds  and  parameters  to  adjust,  and  each 
parameter  should  have  a  distinctly  different  effect  on  the  output.  Here,  we  will 
concentrate  on  the  three  most  important  parameters  discussed  earlier  in  Section  2.2.3, 
namely: 

1)  The  minimum  coherence  threshold,  K ,  which  is  used  to  screen  less  reliable  input 
data  from  the  inversion  process; 

2)  The  Lagrange  multiplier,  Av,  which  controls  the  degree  of  vertical  spectral 
regularization  of  the  input  data;  and, 

3)  The  Lagrange  multiplier,  Ah,  which  controls  the  degree  of  horizontal  spatial 
regularization  of  the  data  filtered  in  Step  2. 

These  three  parameters  are  the  most  important  in  defining  the  operating  points  of  the 
inversion.  Significant  attention  will  be  paid  to  how  the  k  threshold  controls  the  solution 
output.  The  Lagrange  multiplier  for  the  spectral  regularization  will  be  examined  with 
respect  to  its  effect  on  the  vertical  smoothing.  The  Lagrange  multiplier  for  the  horizontal 
regularization  will  be  given  the  least  attention  since  it  amounts  to  a  conceptually  simple 
adaptation  of  the  widely  used  and  well  understood  Bostick  filtering  method. 
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Figure  38.  Complex  Resistivity  Magnitude  Image  using 
Parameters  k  =  .89,  A,v  =  .4,  A,h  =  190. 
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Figure  39.  Complex  Resistivity  Phase  Image  using 
Parameters  K  =  .89,  A,v  =  .4,  A,h  =  190. 
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In  all  cases  where  we  compare  the  effects  of  different  parameter  settings,  we  will  plot  the 
magnitude  and  phase  of  the  complex  resistivity  pc  defined  in  Section  2.2.3,  because  the 
effects  of  the  different  parameter  settings  are  most  clearly  visible  in  pc .  It  is  important  to 
realize  that  the  magnitude  of  pc  alone  is  not  sufficient  to  determine  whether  a 
perturbation  is  conductive  or  resistive  in  nature;  it  is  also  necessary  to  examine  its  phase. 
Once  the  analysis  is  complete,  we  can  plot  an  estimate  of  the  actual  ground  resistivity  as  a 
function  of  depth  taking  into  account  both  the  magnitude  and  phase  of  pc . 

We  turn  first  to  the  complex  resistivity  images  shown  in  Figures  38  and  39.  The 
magnitude  image  shows  a  distinct  anomaly  starting  at  a  frequency  of  about  12  kHz 
located  beneath  stations  8  and  9.  The  peak  of  the  anomaly  occurs  at  a  frequency  of  9  kHz, 
corresponding  to  a  depth  of  about  36  meters  using  the  simple  method  outlined  below  to 
estimate  the  depth.  These  parameters  compare  well  with  the  known  location  of  the  tunnel, 
which  crosses  the  survey  line  between  stations  8  and  9  at  a  distance  of  0.8  meters  from 
station  9,  and  with  the  known  tunnel  depth  of  29  meters. 

The  image  shown  in  Figure  38  was  computed  using  a  minimum  coherence  threshold  of 
k  =.89,  a  spectral  regularization  parameter  of  Av=.4  and  a  Lagrange  multiplier  for  the 
horizontal  of  \=  190.  This  image  is  the  current  best  result,  obtained  by  processing 
remote  reference  impedance  estimates  using  the  APTI  spectral  regularization  technique, 
followed  by  a  modified  Bostick  inversion.  A  true  resistivity  image  corresponding  to  the 
complex  resistivity  images  shown  here  will  be  given  in  Figure  52  at  the  end  of  this 
section. 

The  results  in  Figures  38-39  were  obtained  by  choosing  an  acceptable  operating  point  (k, 
Av  ,  Ah )  subjectively,  not  by  having  objectively  arrived  at  this  choice  as  a  result  of  some 
initial  optimization,  or  by  having  had  taken  advantage  of  some  prior  constraints  resulting 
from  some  trustworthy  prior  knowledge.  Nevertheless,  having  chosen  K ,  Av,  Ah  by 
whatever  means,  we  may  then  discuss  an  optimal  solution  p*  that  minimizes  the 
objective  J[x\Av  ,Ah](pc) ,  determined  by  these  parameters. 
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Many  inversion  algorithms,  including  widely  used  magnetotelluric  inversion  methods, 
have  similar  subjectively  chosen  parameters.  In  this  case,  we  are  attempting  to  estimate  a 
temporally  constant  state  in  which  the  forward  observer  has  a  nontrivial  kernel  that  is  also 
rendered  uncertain  due  to  noise.  In  such  cases,  the  use  of  asymptotic  methods  is 
appropriate.  There  are  many  mathematical  techniques  available  to  address  this  important 
part  of  the  problem,  and  examination  of  the  alternatives  is  an  important  part  of  APTI’s 
algorithm  research.  For  the  moment,  we  will  focus  on  the  more  immediate  aspect  of  the 
problem,  namely  the  design  and  evaluation  of  objectives  such  as  J[/c,Av,/lh](-),  rather 
than  the  feedback  aspect  related  to  J(pc). 

Returning  to  Figure  38,  we  also  observe  a  less  pronounced  anomaly  beneath  stations  2 
and  3  to  the  left  of  the  Golden  Zone  anomaly.  It  is  important  to  understand  that  it  is 
proper  to  refer  to  both  as  resistivity  anomalies  but  we  have  yet  to  discuss  their  phase, 
which  will  determine  whether  they  are  more  or  less  resistive  than  the  background.  Figure 
39  shows  the  phase  corresponding  to  the  magnitude  image  of  Figure  38.  Here,  we  see  that 
the  two  anomalies  yield  a  phase  response  of  opposite  sign.  For  the  case  of  a  void,  one 
expects  that  a  good  estimate  of  the  complex  resistivity  pc  should  have  negative  phase  at 
frequencies  corresponding  to  skin  depths  near  its  actual  depth,  while  a  conductive 
anomaly  should  have  positive  phase.  This  can  be  seen  intuitively  in  two  limiting  cases. 


In  a  homogeneous  medium  the  skin  depth  is  related  to  the  wavelength  according  to 
A  =  2 7tS,  so  that  wavelength  and  target  depth  can  be  related  using  the  approximation 


A 

4 


n  „  nJl  „  ^  „ 

-S  = - 2D  ~ 2D, 

2  4 


where  we  have  used  the  Bostick  estimate  of  the  target  depth,  D=S/-J2.  This  means  that 
a  surface  field  excitation  has  a  90°  phase  lead  over  the  surface  response,  if  there  is  any, 
due  to  a  target  at  depth  D.  In  the  case  of  a  target  satisfying  the  relatively  deep 
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approximation,  namely  that  its  depth  and  diameter  are  in  the  relation  D»  d  ,  the  surface 
response  for  a  TM  mode  excitation  is  almost  entirely  due  to  static  charge  collection  at  the 
target  surface,  and  is  well  modeled  as  an  equivalent  two-dimensional  dipole.  In  the  case 
of  a  void,  this  charge  collection  leads  to  a  secondary  electric  field  that  is  in-phase  with  the 
primary  electric  field  at  depth.  As  a  result,  the  secondary  electric  field  at  the  surface  lags 
the  primary  magnetic  field  in  phase  by  45°  and  therefore  produces  a  pc  that  can  be 

approximated  as 


p0co 


If  E0+E1 


\2 


Hn 


so  that  its  phase  satisfies 


tan  Zpc  ~ 


from  which  it  follows  that  tan  Zpc  <  0  and  in  the  limit  of  small  response  we  have 


tan  Zpe  = 


In  other  words,  the  complex  resistivity  will  exhibit  at  most  a  90°  phase  lag,  similar  to  a 
capacitor  connected  in  parallel,  as  a  result  of  the  in-phase  charge  collection  on  the  void. 
In  the  case  of  a  conductive  body,  the  secondary  electric  field  is  180°  out-of-phase  with 
the  primary  at  the  surface  of  the  target  and  this  results  in 
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so  that  the  complex  resistivity  will  exhibit  an  inductive  phase  lead  of  at  most  90°.  Of 
course,  the  above  analysis  employs  a  single  lumped  circuit  argument  that  only  grossly 
approximates  what  is  in  truth  a  three-dimensional  ground.  For  completeness,  we  mention 
that  a  similar  lumped  circuit  analysis  can  be  applied  when  the  target  satisfies  a  relatively 
shallow  approximation.  Under  this  assumption,  two-layer  half  space  models  apply  — 
relatively  shallow  voids  and  conductors  corresponding  to  more  conductive  and  less 
conductive  top  layers,  respectively  —  and  can  be  used  to  show  that  the  same  phase 
characteristics  are  obtained. 

In  general,  the  lead  or  lag  characteristics  of  the  complex  resistivity  response  due  to  the 
extremes  of  conductive  targets  or  resistive  voids  can  be  expected  to  be  less  than  the  ±90° 
phase  differences  of  ideal  circuit  elements.  Nevertheless,  significant  and  consistent 
departures  of  the  phase  of  the  complex  resistivity  away  from  zero  indicate  the  conductive 
or  resistive  nature  (relative  to  the  background)  of  anomalies  that  are  observed. 

Returning  to  Figure  38,  we  see  that  the  resistivity  anomaly  on  the  right,  which  correlates 
well  with  the  known  horizontal  position  and  depth  of  the  Golden  Zone  tunnel,  exhibits  a 
negative  phase  anomaly  as  seen  in  Figure  39.  This  is  consistent,  in  light  of  the  above 
analysis,  with  the  presence  of  a  void.  On  the  other  hand,  Figure  39  appears  to  indicate 
that  the  anomaly  to  the  left,  beneath  stations  2  and  3,  is  conductive  in  nature. 

We  turn  next  to  the  examination  of  the  effect  of  changing  the  minimum  coherence 
threshold  k  ,  as  discussed  in  Section  2.2.3,  for  fixed  values  of  and  \ .  A  statistic  that 
is  useful  when  adjusting  the  coherence  threshold  is  the  percentage  of  the  input  data  that 
exceed  the  threshold.  As  a  first  example,  the  coherence  threshold  of  K=. 89  that  was  used 
in  computing  Figure  38  resulted  in  2.54%  of  the  input  data  having  a  coherence  greater 
than  .89;  as  a  shorthand  we  will  write  this  as  %(.89)=2.54.  Figures  40  through  44  give  the 


90 


results  of  increasing  K  while  keeping  the  other  parameters  fixed  at  Av=.4  and  A^  =  190, 
as  were  used  in  the  first  two  figures.  For  brevity,  we  will  show  only  the  magnitude 
images,  although  complete  interpretation  of  the  data  would  require  the  phase  as  well. 

Figure  40  is  important  in  that  it  has  its  minimum  coherence  threshold  set  to  zero,  so  that 
%(0)=100;  in  other  words,  all  the  input  data  at  every  frequency  was  included  in  the 
inversion.  Although  remote  reference  processing  was  employed,  so  that  the  input 
impedance  estimates  were  unbiased,  the  peaks  that  can  be  seen  in  the  resistivity  image 
appear  to  be  too  high  in  magnitude,  whether  they  are  due  to  an  underlying  ground 
anomaly  or  not.  This  suggests  some  future  attention  be  paid  to  robust  remote  reference 
impedance  estimation. 

Figures  41  through  44  present  the  remainder  of  the  parameter  scan  through  the  five  values 
of  k  -  80,  85,  87,  and  95,  respectively,  which  may  be  compared  to  the  other  figures;  the 
associated  percentages  here  are  given  by  %(85)=3.9,  %(87)=3.0,  %(95)=1.8.  We  observe 
that  the  anomaly  associated  with  the  tunnel  is  obscured  by  the  inclusion  of  the  relatively 
noisy  data  that  passes  the  lower  coherence  thresholds,  but  appears  strongly  when  a  high 
threshold  is  set.  This  effect  is  due  to  the  abrupt  nature  of  the  thresholding  process. 

It  should  be  mentioned  that  many  MT  inversion  algorithms  that  employ  remote  reference 
processing  use  no  K  thresholding,  since  the  cross-spectra  from  remote  reference 
processing  are  unbiased  estimates  of  the  true  spectra  regardless  of  the  coherence. 
However,  the  coherence  is  a  measure  of  the  degree  to  which  the  measured  data  are 
consistent  with  a  distant  plane  wave,  and  so  the  coherence  may  be  used  with  remote 
reference  data  as  an  indication  of  the  quality  of  the  data. 
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Figure  40.  Complex  Resistivity  Magnitude  Image  using 

Parameters  k  =  0,  E,  =  A,  lh  =  190;  %(0)  =  100. 
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Figure  41 .  Complex  Resistivity  Magnitude  Image  using 

Parameters  k  =  .80,  Av  =  .4,  Ah  =  190;  %(80)  =  8.2. 
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Figure  42.  Complex  Resistivity  Magnitude  Image  using 

Parameters  k  =  .85,  Xv  =  .4,  A,h  =  190;  %(85)  =  3.96. 
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Figure  43.  Complex  Resistivity  Magnitude  Image  using 

Parameters  k  =  .87,  Xv  =  .4,  A,h  =  190;  %(87)  =  3.01. 
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Figure  44.  Complex  Resistivity  Magnitude  Image  using 

Parameters  k  =  .95,  Av  =  .4,  =  190;  %(95 )  =  1.80. 
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We  have  seen  that  the  abrupt  inclusion  or  exclusion  of  data  by  coherence  thresholding 
makes  the  images  very  sensitive  to  the  value  of  K .  Ideally,  a  continuous  weighting  would 
be  employed  to  reduce  the  effects  of  poor  quality  data  without  excluding  it  entirely.  This 
would  reduce  the  sensitivity  of  the  image  to  settings  of  the  weighting  parameter.  We  have 
employed  the  coherence  here  as  an  easily  computed  indicator  of  data  quality,  but  it  is  also 
possible  to  derive  more  accurate  quality  estimates  using  the  electric  and  magnetic  field 
signal-to-noise  ratios  that  can  be  computed  using  a  remote  reference.  Implementation  of 
this  enhancement  to  the  processing  procedure  is  planned  for  the  next  year. 

An  additional  effect  of  changes  in  the  threshold  is  that  the  average  conductivity  decreases 
as  the  K  threshold  increases.  This  is  because  as  the  percentage  of  points  that  are  used  in 
the  inversion  falls,  more  weight  is  implicitly  given  to  the  prior  constraint. 

Directly  modifying  results  in  a  similar  bias.  In  this  sense,  the  complex  resistivity 
images  shown  are  more  properly  considered  as  indicators  of  significant  departures  from  a 
one-dimensional  stratified  earth  (that  is,  indicators  of  the  presence  of  two-  and  three- 
dimensional  objects),  rather  than  estimates  of  the  background  resistivity.  Conversion  of 
the  complex  resistivity  to  an  accurate  background  resistivity  estimate  and  production  of  a 
depth  estimate  consistent  with  that  resistivity  is  the  final  step  in  producing  an  image,  as 
will  be  described  at  the  end  of  this  section. 
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Figure  45.  Complex  Resistivity  Magnitude  Image  using 
Parameters  k  =  .89,  Xv  =  .05,  A,h  =  190. 
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Figure  46.  Complex  Resistivity  Magnitude  Image  using 
Parameters  k  =  .89,  Xv  =  .20,  A,h  =  190. 
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Figure  47.  Complex  Resistivity  Magnitude  Image  using 
Parameters  k  =  .89,  A,v  =  .80,  A,h  =  190. 
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Figure  48.  Complex  Resistivity  Magnitude  Image  using 
Parameters  k  =  .89,  =  6.4,  =  190. 
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Next  we  turn  to  examination  of  p*(. 89,  Av,  190)  as  Av  is  varied.  Figures  45  -  48  give  the 
magnitude  and  phase  images  for  the  values  /tv=.05,  0.2,  0.8  and  6.4.  We  see  that  as  Av 
increases  in  value,  the  estimate  displays  less  resolution  in  the  vertical  direction,  which  is 
a  reflection  of  the  fact  that  the  Jv[/c,Av](/?c)  objective  determines  an  equivalent  vertical 
filter  width  that  becomes  broader  at  a  given  frequency  as  Av  increases.  On  the  other 
hand,  if  Av  is  set  too  small  (compare  Av=.05  and  0.2),  then  there  may  be  insufficient 
support  for  the  vertical  filter  at  a  given  frequency  to  achieve  constructive  superposition. 

Figures  49  through  51  show  p*(. 89,  .4,  Ah)  as  Ah  is  varied  to  Ah  =120,  220  and  320. 
The  effect  of  this  straightforward  adaptation  of  the  Bostick  inverse  is  easy  to  appreciate, 
in  that  it  simply  leads  to  a  loss  of  horizontal  resolution  as  Ah  is  raised,  and  for  a  fixed 
Ah  ,  a  resolution  loss  that  increases  with  depth. 
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Figure  49.  Complex  Resistivity  Magnitude  Image  using 
Parameters  k  =  .89,  A,v  =  .4,  A,h  =  120. 
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Figure  50.  Complex  Resistivity  Magnitude  Image  using 
Parameters  k  =  .89,  =  .4,  A,h  =  220. 
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Figure  51 .  Complex  Resistivity  Magnitude  Image  using 
Parameters  k  =  .89,  Xv  =  .4,  =  320. 
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4.3.3.  Final  Resistivity  Image 


We  next  turn  the  discussion  to  the  estimation  of  a  real-valued  resistivity  image,  as 
opposed  to  the  complex  resistivity  images  that  has  been  the  focus  until  now.  As  we  have 
outlined  earlier,  the  TM  response  of  a  two-dimensional  void  satisfying  the  relatively  deep 
approximation,  D»d ,  is  well  modeled  as  an  equivalent  quasi-static  dipole  in-phase 
with  the  current  excitation  at  depth.  For  this  case,  in  the  limit  defining  the  electrically 
shallow  approximation,  i.e.,  D«S(f),  we  see  that  the  current  due  to  the  equivalent 
dipole  response  at  the  surface  will  be  in-phase  with  the  surface  current  due  to  the  primary 
field.  In  this  case,  a  void  will  result  in  a  distinct  rise  in  the  magnitude  of  the  complex 
resistivity  while  a  metallic  cylinder  will  result  in  its  decrease,  as  one  might  expect. 

However,  in  the  discussion  above,  we  have  shown  that  the  well  known  Bostick  apparent 
depth  criterion,  which  is  valid  for  a  half-space,  defines  a  frequency  at  which  the  quarter 
wavelength  approximates  the  round-trip  distance  to  the  target.  At  this  frequency,  for 
which  we  have  only  D  ~  S(f) ,  the  equivalent  quasi-static  dipole  will  yield  a  surface 
complex  resistivity  that  has  an  increase  in  magnitude  for  either  a  resistive  void  or  a 
conductive  cylinder.  The  two  cases  will,  however,  have  opposite  complex  resistivity 
phase  responses. 

From  the  complex  resistivity  magnitude  and  phase,  we  can  compute  a  final  output  image 
that  includes  the  information  from  both  quantities  and  constitutes  our  best  estimate  of  the 
actual  (non-complex)  subsurface  resistivity.  In  addition,  we  can  compute  depth  estimates 
for  the  frequencies  at  which  these  real  resistivities  are  computed,  resulting  in  an  image 
showing  the  estimated  subsurface  resistivity  as  a  function  of  depth.  This  image  is  shown 
in  Figure  52. 
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Figure  52.  Final  Resistivity  Image  of  Golden  Zone 
Line  1  Data. 
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The  resistivity  estimate  in  Figure  52  is  derived  from  the  complex  resistivity  in  the 
following  manner.  As  discussed  above,  imaging  in  the  TM  case  involves  a  perturbation 
solution  for  the  problem,  involving  relatively  deep  loci  of  quasi-static  charge.  This 
perturbation  solution  may  be  easily  recast  in  terms  of  a  positive,  real-valued  resistivity  if 
we  superpose  the  magnitude  and  phase  anomalies  in  the  complex  resistivity  pc  using  the 


expression 


( 


V 


\ 

sin  /-pc 

J 


where  pb  is  an  estimate  of  the  background  resistivity  obtained,  for  example,  using 
classic,  low-resolution  resistivity  estimation  techniques. 

The  final  image  is  consistent  with  most  known  information  on  the  characteristics  of  the 
subsurface  at  survey  line  1.  A  partially  resolved  surface  conducting  layer  is  present,  as  is 
common,  although  it  is  not  clear  that  the  surface  layer  at  survey  line  1  extended  to  the  15 
meter  depth  shown.  The  site  had  been  terraced  by  earthmoving  equipment,  so  it  is 
possible  that  the  disturbed  layer  extended  to  this  depth.  It  is  more  likely,  however,  that  the 
conducting  layer  is  not  completely  resolved,  and  that  a  thinner  layer  of  higher 
conductivity  has  been  imaged  as  a  thick  layer  because  the  frequency  data  above  50  kHz 
needed  to  resolve  the  thin  layer  is  not  present  in  the  input. 

The  location  and  depth  of  the  resistive  anomaly  shown  beneath  stations  8  and  9  are 
consistent  with  the  known  location  of  the  tunnel.  The  horizontal  location  of  the  anomaly 
matches  the  tunnel  location  to  within  2-3  meters,  approximately  the  limit  of  horizontal 
accuracy  for  a  station  spacing  of  5  meters.  The  estimated  depth  of  34  meters  is  larger  than 
the  known  tunnel  depth  of  29  meters.  Accurate  estimation  of  depths  when  a  partially 
resolved  surface  layer  is  present  is  very  difficult,  and  errors  of  20  percent  are  not 
unexpected.  (A  more  complete  discussion  of  the  depth  estimation  procedure  is  given 
below.) 
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The  conductive  anomaly  seen  beneath  stations  2  and  3  is  of  unknown  origin.  It  appears  at 
approximately  the  same  depth  as  the  resistive  anomaly  believed  to  be  caused  by  the 
tunnel,  and  is  of  relatively  small  size.  This  suggests  that  it  might  be  due  to  a  man-made 
object,  but  it  is  unknown  whether  tunneling  occurred  at  this  level  before  the  current  325 
foot  level  tunnel  was  constructed  in  1983. 

The  background  resistivity  estimated  for  the  region  beneath  survey  line  1  is  80-100  Ohm- 
meters.  Although  there  are  no  measurements  of  this  exact  location  for  comparison,  this  is 
somewhat  lower  than  the  200-300  Ohm-meters  computed  for  the  general  area  from  earlier 
airborne  loop-loop  survey  measurements  at  9  kHz.  Since  both  the  airborne  survey  and  this 
experiment  measured  resistivities  in  approximately  the  same  depth  range,  better 
agreement  between  the  two  estimates  would  have  been  expected.  The  resolution  of  the 
airborne  survey  data  were  rather  coarse,  however,  and  it  is  difficult  to  determine  whether 
local  variations  in  conductivity  that  were  not  resolved  by  the  airborne  data  are  sufficiently 
large  to  account  for  this  difference.  Processing  of  the  data  from  survey  lines  2  and  3 
should  provide  some  additional  information  to  address  this  question. 

4.3.4.  Depth  Estimation 

We  will  now  describe  the  depth  estimation  procedure  for  the  final  resistivity  image  in 
more  detail.  It  should  be  mentioned  that  depth  estimation  is  often  treated  as  a  decoupled 
problem,  i.e.,  given  a  resistivity  estimate,  employ  it  directly  to  estimate  the  depth.  In  the 
general  case,  the  two  should  be  estimated  simultaneously.  This  is  especially  true  since, 
given  its  nature,  one  knows  a  priori  that  any  depth  estimate  should  be  monotonic. 

Direct  employment  of  this  constraint,  even  in  a  decoupled  approach,  is  somewhat  difficult 
since  its  description  leads  to  nonlinearities  in  the  design  of  a  constraint  objective.  In 
general,  the  coupled  depth-resistivity  problem  deserves  and  requires  more  attention, 
however,  at  present  we  employ  a  simplified  approach  that  can  be  trusted,  and  has  been 
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shown,  to  provide  reasonable  estimates  in  relatively  simple  geologies.  As  a  compromise, 
our  approach  addresses  the  monotonicity  aspect  of  the  depth  estimate  in  a  manner  that  is 
coupled  with  a  resistivity  estimate,  but  avoids  nonlinearities. 

Put  simply,  the  depth  averaging  algorithm  begins  by  taking  the  average  of  the  remote 
reference  impedance  at  all  measured  frequencies  using  a  very  low  coherence  threshold. 
We  used  k=.1  and  in  addition  deleted  all  data  for  which  the  frequencies  displayed  near¬ 
field  effects.  This  impedance  average  was  then  processed  using  the  well-known  Cagniard 
method  which  produces  a  stable  resistivity  estimate  when  the  ground  is  close  to  one¬ 
dimensional.  To  address  possible  errors  in  the  estimate  due  to  this  one-dimensional 
assumption,  the  Cagniard  output  is  next  processed  to  obtain  a  marginal  using  a  classical 
median  filter  at  each  frequency.  This  one-dimensional  marginal  delivers  what  may  be 
considered  a  “background”  resistivity  for  the  region  in  question. 

Next,  the  resistivity  marginal  is  subjected  to  an  iterative  smoothing  process  to  enforce  the 
monotonicity  constraint.  As  an  initial  condition,  the  resistivity  marginal  is  used  to  obtain 
a  point-wise  estimate  of  the  skin  depth  using  the  Cagniard  method.  This  initial  estimate  is 
not  in  general  monotonic.  Then,  the  resistivity  marginal  is  successively  filtered  using  a 
low-pass  FIR  filter  having  a  relatively  wide  bandwidth,  until  the  skin  depth  computed 
using  the  Cagniard  method  satisfies  a  monotonicity  constraint  to  within  a  desired 
tolerance.  The  skin  depth  estimate  which  results  is  used  to  assign  depths  to  the 
resistivities  estimated  above,  resulting  in  the  depth  scale  shown  in  the  final  image. 
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5.  Conclusions  and  Future  Plans 


The  results  described  in  the  preceding  sections  include  several  important 
accomplishments,  including  design  and  integration  of  the  hardware  and  software  needed 
for  collecting  field  data,  detailed  electromagnetic  simulations  that  confirm  several 
important  effects  seen  in  experimental  data,  and  a  successful  field  experiment  that 
resulted  in  imaging  of  the  Golden  Zone  mine  tunnel.  Major  areas  of  emphasis  for  the  next 
year  of  this  effort  include: 

(1)  Elimination  of  the  backlog  of  unprocessed  experimental  data.  This  backlog 
results  from  both  the  data  transfer  difficulties  described  in  Section  4. 1.2.5,  and  from  the 
lack  of  a  packaged  inversion  algorithm  that  could  produce  images  without  input  from  an 
operator.  Both  of  these  problems  are  in  the  process  of  being  resolved,  and  it  is  expected 
that  the  backlog  will  be  eliminated  within  the  next  three  months. 

(2)  Improvement  of  the  speed  with  which  new  data  are  processed.  The  same 
modifications  implemented  to  remove  the  backlog  of  unprocessed  data  will  also  allow 
much  faster  processing  of  new  data,  so  that  the  results  of  experiments  are  available  within 
weeks  after  the  completion  of  the  experiment. 

(3)  More  extensive  use  of  simulation  to  study  imaging  performance.  Current 
simulations  are  not  exact  replacements  for  experimental  data,  principally  because  they 
contain  far  fewer  frequencies:  typical  simulations  have  data  at  tens  or  hundreds  of 
frequencies,  rather  than  the  thousands  in  an  experimental  data  set.  Production  of 
simulated  data  sets  that  match  actual  experiments  in  frequency  resolution  will  allow 
simulated  data  to  substitute  for  experimental  results  in  imaging  studies. 

(4)  Completion  of  detailed  comparisons  between  simulated  data  and  experimental 
data,  and  between  APTI’s  inversion  algorithm  and  conventional  magnetotelluric 
algorithms.  The  availability  of  high  frequency  resolution  simulations  will  also  make  these 
comparisons  possible,  allowing  us  to  verify  conclusions  drawn  from  experimental  data 
and  examine  the  improvements  in  imaging  performance  possible  using  APTI’s  imaging 
techniques. 
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(5)  Conduct  field  experiments  on  a  wider  variety  of  targets,  including  conductive 
targets.  The  experiments  conducted  to  date  have  concentrated  on  simple  two-dimensional 
resistive  tunnels.  Many  other  types  of  structures  are  also  of  military  interest,  and  the  other 
simple  two-dimensional  target,  a  conductive  object,  will  also  be  examined  in  experiments 
and  computations.  Since  the  field  perturbations  at  the  surface  due  to  this  type  of  target  are 
larger  than  those  due  to  a  resistive  object,  it  is  expected  that  better  performance  can  be 
achieved  on  conductive  targets.  « 
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